Metodos Numericos Examen
Metodos Numericos Examen
Metodos Numericos Examen
FUNDAMENTOS MATEMATICOS
DE LA INGENIERIA II
ETSIT
Angel
Duran.
Departamento de Matematica Aplicada.
Universidad de Valladolid
14 de enero de 2011
Contenidos
1. Preliminares
1.1. Introduccion . . . . . . . . . . . . . . . . . . . . . .
1.2. El concepto de algoritmo . . . . . . . . . . . . . . .
1.3. Aritmetica del ordenador . . . . . . . . . . . . . .
1.3.1. Representacion binaria . . . . . . . . . . . .
1.3.2. Representacion en el ordenador . . . . . . .
1.4. Analisis del error . . . . . . . . . . . . . . . . . . .
1.4.1. Redondeo por aproximacion y redondeo por
1.4.2. Tipos de errores y propagacion del error . .
1.5. Ejercicios . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
4
4
5
6
6
9
11
11
13
18
21
21
23
23
26
28
29
31
36
40
40
44
46
49
57
60
62
62
68
72
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
eliminacion (o truncamiento)
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
3. Interpolaci
on polin
omica
74
3.1. Introduccion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.2. Interpolacion de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
3.3.
3.4.
3.5.
3.6.
4. M
etodos iterativos
4.1. Introduccion . . . . . . . . . . . . . . . . . . .
4.2. Metodos iterativos para sistemas lineales . . .
4.2.1. Tecnicas iterativas clasicas . . . . . .
4.2.2. Metodo de Jacobi . . . . . . . . . . .
4.2.3. Metodo de Gauss-Seidel . . . . . . . .
4.2.4. Metodo SOR . . . . . . . . . . . . . .
4.2.5. Analisis de convergencia . . . . . . . .
4.3. Metodos iterativos para ecuaciones no lineales
4.3.1. Metodo de biseccion . . . . . . . . . .
4.3.2. Metodo de la secante . . . . . . . . . .
4.3.3. Metodo de punto fijo . . . . . . . . . .
4.3.4. El metodo de Newton . . . . . . . . .
4.3.5. Aspectos de convergencia . . . . . . .
4.4. Metodos iterativos para sistemas no lineales .
4.4.1. Iteracion de punto fijo para sistemas .
4.4.2. El metodo de Newton para sistemas .
4.5. Ejercicios . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
76
81
82
84
85
85
86
89
93
97
100
101
106
109
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
111
111
113
113
113
115
117
118
124
124
126
128
133
134
136
136
136
137
.
.
.
.
.
.
.
.
.
.
.
.
143
. 143
. 144
. 145
. 147
. 148
. 149
. 151
. 152
. 153
. 155
. 158
. 158
5. Cuadratura y derivaci
on num
erica
5.1. Introduccion . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2. Metodos elementales de construccion de reglas de cuadratura
5.2.1. Metodo interpolatorio . . . . . . . . . . . . . . . . . .
5.2.2. Metodo directo . . . . . . . . . . . . . . . . . . . . . .
5.2.3. Metodo del desarrollo de Taylor . . . . . . . . . . . . .
5.3. Analisis del error . . . . . . . . . . . . . . . . . . . . . . . . .
5.4. Reglas compuestas . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1. Analisis del error en las reglas compuestas . . . . . . .
5.4.2. Estudio comparativo . . . . . . . . . . . . . . . . . . .
5.5. Formulas de aproximacion a la derivada . . . . . . . . . . . .
5.5.1. Aproximacion a derivadas de orden superior . . . . . .
5.5.2. Error en la derivacion numerica . . . . . . . . . . . . .
2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
164
164
166
167
172
178
188
188
189
191
194
195
200
203
204
Tema 1
Preliminares
1. Introduccion.
2. Concepto de algoritmo.
3. Aritmetica del ordenador.
4. Analisis del error.
5. Ejercicios.
1.1.
Introducci
on
1.2.
El concepto de algoritmo
k = N 1, . . . , 0,
entonces b0 = P (c)
El algoritmo tiene ademas relacion con el proceso de division de polinomios (veanse los ejercicios).
Es u
til familiarizarse con una escritura de los algoritmos independiente del lenguaje de programacion utilizado, que llamaremos pseudocodigo. Mostrara claramente los pasos de que consta
el algoritmo, de manera que su traslacion a cualquier lenguaje (FORTRAN, C, MATLAB, etc)
5
aN
c
bN
aN 1
cbN
bN 1
ak
cbk+1
bk
a2
cb3
b2
a1
cb2
b1
a0
cb1
b0 = P (c)
sea directa. El pseudocodigo debe incluir ademas aspectos propiamente computacionales del algoritmo. Por ejemplo, un pseudocodigo para el algoritmo de Horner se escribira:
Pseudoc
odigo del algoritmo de Horner
Dados a0 , . . . , aN y c,
Desde k = N 1 hasta 0
ak ak + cak+1
Puede observarse el ahorro computacional al almacenar cada parentesis en el propio coeficiente
ak , de manera que el valor del polinomio en c esta finalmente contenido en a0 .
Puede haber diferentes algoritmos que resuelvan un mismo problema. Y varios pueden ser
los criterios para sus discriminacion [7]. El mas habitual es el que selecciona aquel metodo que,
permitiendo que se generen unos errores limitados a priori, realice menos operaciones y, por tanto,
el ordenador tarde menos tiempo en su ejecucion. Tanto la generacion de errores como el coste
computacional u operativo constituyen lo que se llama el an
alisis de eficiencia del algoritmo.
Una medida del coste computacional suele ser el n
umero de operaciones. De estas, suelen
tener mas importancia las multiplicaciones y divisiones, pues estas operaciones son mas costosas
en tiempo para el ordenador que las sumas y restas.
Ejemplo 1.2. Por ejemplo, el algoritmo de Horner puede compararse con el procedimiento
elemental de evaluar P (x) en x = c a traves de las potencias de c, c2 = c c, c3 = c2 c, . . . Este
u
ltimo requiere N 1 multiplicaciones para las potencias de c, otras N multiplicaciones para
los productos aj cj , j = 1, . . . , N y N sumas para la combinacion lineal. Esto da un total de
2N 1 multiplicaciones y N sumas. Por otro lado, el algoritmo de Horner (se puede contar en el
pseudocodigo) requiere N multiplicaciones y N sumas. De este modo, el algoritmo de Horner es
mas eficiente, algo que se notara mas si el grado del polinomio es muy alto o el polinomio ha de
evaluarse en muchos puntos.
1.3.
Aritm
etica del ordenador
Hay que tener en cuenta que, con independencia del algoritmo, su implementacion en el
ordenador ha de considerar el hecho de que este utiliza aritmetica finita, con la consiguiente
introduccion de errores. Debemos entonces comentar algo sobre la representacion de los n
umeros
en el ordenador.
1.3.1.
Representaci
on binaria
natural N , la representacion de un n
umero racional T en base N es como sigue: si ai es el dgito
situado en la i-esima posicion
T = (ap ap1 a2 a1 a0 a1 a2 aq )N = ap N p + ap1 N p1 + + aq N q .
Los dgitos ap , . . . , aq varan entre los N n
umeros diferentes {0, 1, . . . , N 1}. Por ejemplo
73,35 = (73,35)10 = 7 101 + 3 100 + 3 101 + 5 102 ,
(10,11)2 = 1 21 + 0 20 + 1 21 + 1 22 .
Decimal
Base 10
0
1
1
3
4
5
6
7
8
9
10
11
12
13
14
15
Binaria
Base 2
0
1
10
11
100
101
110
111
1000
1001
1010
1011
1100
1101
1110
1111
Octal
Base 8
0
1
2
3
4
5
6
7
10
11
12
13
14
15
16
17
Hexadecimal
Base 16
0
1
2
3
4
5
6
7
8
9
A
B
C
D
E
F
= b1 N 1 + + bs+1 N s+1 + bs N s
= N 1 (b1 + N 1 (b2 + r + N 1 (bs+1 + N 1 bs ) . . .),
Fk = f rac(N Fk1 ),
k 1.
bj N j ,
bj {0, 1, . . . , N 1}.
j=1
Pseudoc
odigo del algoritmo 2
Dado T = E + F y N
Poner Q0 = E, j = 0
while Q0 6= 0
Q1 = ent(Q0 /N )
bj = Q0 N Q1
j =j+1
Q0 Q1
end while
Ejemplo 1.3. Para representar el n
umero T = 33 en el sistema binario, se tiene el algoritmo
33 = 2 16 + 1,
b0 = 1,
16 = 2 8 + 0,
b1 = 0,
8 = 2 4 + 0,
b2 = 0,
4 = 2 2 + 0,
b3 = 0,
2 = 2 1 + 0,
b4 = 0,
1 = 2 0 + 1,
b5 = 1.
b1 = 1,
F1 = 0,4,
2 F1 = 0,8,
b2 = 0,
F2 = 0,8,
2 F2 = 1,6,
b3 = 1,
F3 = 0,6,
2 F3 = 1,2,
b4 = 1,
F4 = 0,2,
2 F4 = 0,4,
b5 = 0,
F5 = 0,4,
2 F5 = 0,8,
b6 = 0,
F6 = 0,8,
2 F6 = 1,6,
b7 = 1,
F7 = 0,6,
1.3.2.
Representaci
on en el ordenador
dS1
dS2
dS1
d2 d1
0 +
=
1
d0
1 2S1
= 2S1 1.
12
De este modo, con S posiciones, se pueden almacenar los enteros entre (2S1 1) y (2S1 1).
Los valores habituales utilizados son S = 16 (con |Nmax | = 32767) y S = 32 (con |Nmax | =
2147483647).
Almacenamiento de n
umeros reales
Los ordenadores almacenan los n
umeros reales en una representacion binaria llamada de
coma flotante normalizada [3]. Esto significa que lo que se almacena no es el n
umero x sino una
aproximacion binaria a x,
x M 2E ,
donde M se llama mantisa y E exponente. Esta representacion servira para cualquier base de numeracion, no solo la binaria. Por ejemplo, el n
umero (891,246)10 se representa como (0,891246)10
3
10 (mantisa 0,891246, exponente 3). El n
umero binario (101,001)2 es (0,101001)2 2(11)2 (mantisa (0,101001)2 , exponente (11)2 ). As, el punto (o coma) decimal se desplaza (flota) hasta que
la parte entera es nula y el primer dgito decimal es no nulo.
Hemos comentado que entonces todos los ordenadores almacenan los n
umeros reales en base
dos y coma flotante. La diferencia entre los ordenadores se encuentra en el n
umero de dgitos
(bits en el sistema binario) para la mantisa M y el exponente E, as como en el orden en que
estos se almacenan.
d1
d2
d(M 1)
dM =
dM
0 +
1
dE1
dE1 =
dE2
0
1
d1
d0
E2 +20 )
10
= (1 2(1M ) )22
E1 1
E2 +d 20 )
0
E2 +20 )
= 22
E1
hay que observar tambien que la mantisa mnima (en base diez y valor absoluto) es 1/2 = 21
(notese que el primer dgito decimal de la mantisa ha de ser no nulo, luego, en base dos, sera (0,1)2 ,
es decir, 1 21 ) y que la mantisa siempre es menor que uno, pues la maxima verifica
21 + 2(M 1) <
X
1
= 1.
2n
n=1
1.4.
An
alisis del error
Ep
.
p
1.4.1.
Sea p un n
umero real, expresado en forma decimal
p = 0.d1 dk+1 10n ,
11
22
= 3,142857142857142857 . . . ,
7
Por ejemplo, p = 0,9951 tiene tres cifras significativas correctas de p = 0,9949, pues
|p p|
= 0,20103 103 .
p
Un fenomeno relacionado con este asunto es el de la cancelacion o perdida de cifras significativas [3]. Por ejemplo, los n
umeros p1 = 1,234567891, p2 = 1,234567892 tienen 10 cifras significativas, mientras que p2 p1 = 0,000000001 solo tiene una. Esto puede producir una reducci
on en
la precision de los calculos que estemos realizando. Un ejemplo clasico [4] es el siguiente:
Ejemplo 1.6. Vamos a comparar los resultados de calcularf (500) y g(500), usando seis
cifras
=
=
= 11,1748.
g(500) =
22,3830 + 22,3607
44,7437
501 + 500
La respuesta g(500) = 11,1748 tiene un error absoluto menor y es la que se obtendra redondeando
por aproximacion la respuesta exacta 11,174755300747198 . . . a seis cifras significativas. La mala
respuesta de f (500) es debida a la perdida de cifras significativas al realizar la resta de las dos
cantidades.
1.4.2.
En [3], los autores dan la siguiente clasificacion de errores que puede contener un n
umero:
1. Errores inherentes: cometidos en los datos al obtenerse de forma experimental.
2. Errores de redondeo: cometidos al almacenar un n
umero real en el ordenador.
3. Errores de truncamiento: cometidos al adaptar un proceso finito (un algoritmo) a uno
infinito, como, por ejemplo, al aproximar series infinitas.
Pasemos ahora a medir la propagacion del error al realizar operaciones. Sean p, q valores exactos
de dos n
umeros y p, q sus correspondientes aproximaciones. Recordemos que Ep , Rp representan,
respectivamente, los errores absoluto y relativo para el n
umero p (lo mismo para q). Vamos a
buscar cotas de error en los resultados obtenidos al sumar, restar, multiplicar y dividir estas
cantidades, as como al evaluar una funcion.
Propagaci
on en la suma. Si s = p+q es el valor exacto de la suma y s = p+ q su aproximaci
on,
entonces
Es = s s = p + q p q = Ep + Eq ,
mientras que
Rs =
Ep + Eq
Es
p
q
=
=
Rp +
Rq .
s
p+q
p+q
p+q
De manera que el error relativo en la suma, con referencia a los errores relativos para los sumandos,
queda amplificado por las cantidades p/(p + q) y q/(p + q), amplificacion que puede ser notable
si se suman dos cantidades opuestas de valor absoluto muy cercano .
13
Propagaci
on en la resta. Si r = p q es el valor exacto de la diferencia y r = p q su
aproximacion, entonces
Er = r r = Ep Eq ,
Er
p
q
Rr =
=
Rp
Rq .
r
pq
pq
El error relativo en la resta, con referencia a los errores relativos para p y q, queda amplificado
por las cantidades p/(p q) y q/(p q). Tambien puede observarse el efecto de amplificaci
on en
el error si se restan dos n
umeros muy parecidos.
Propagaci
on en el producto. Si x = pq es el valor exacto del producto y x
= pq su aproximacion, entonces
Ex = pq pq = pq (p Ep )(q Eq )
Rx
donde en las expresiones anteriores hemos utilizado que la cantidad Ep Eq es despreciable frente
a las otras. De este modo, el error relativo en el producto es muy cercano a la suma de los errores
relativos de los factores.
Propagaci
on en la divisi
on. Si x = p/q es el cociente exacto y x
= p/
q el aproximado, entonces
Ex =
=
p p p p Ep
=
q q q
q Eq
qEp pEq
qEp pEq
.
q(q Eq )
q2
qEp pEq
Rp Rq
Ex
=
=
Rp R q .
x
p(q Eq )
1 Rq
Propagaci
on en la evaluaci
on de una funci
on. sea z = f (p) la evaluacion exacta de una
funcion f (suficientemente derivable) en el n
umero p y z = f (
p) el valor de f en su aproximaci
on.
Entonces
Ez = f (p) f (p Ep ) = f 0 (p)Ep f 00 (p)
Ep2
+ f 0 (p)Ep ,
2
Ep2
Ez
1
=
(f 0 (p)Ep f 00 (p)
+ )
z
f (p)
2
f 0 (p)
p
Rp .
f (p)
Rz =
Los errores sufren una mayor amplificacion si la derivada de f en p es muy alta o si p esta cerca
de un cero de f .
14
Ejemplo 1.7. Vamos a ilustrar lo expresado anteriormente con uno de los ejemplos proporcionados por [3]. Se desea realizar el calculo
z = a(b + c),
(1.1)
(1.2)
Queremos determinar con cual de las dos expresiones anteriores se obtiene una cota del error
menor. Denotaremos por R los errores relativos de las operaciones para la primera expresi
on y
los errores relativos para la segunda.
por R
Los errores en cada una de las operaciones parciales, siguiendo lo que hemos deducido sobre
la propagacion en las operaciones, son, para el primer caso
b
c
rb +
rc + r1 ,
b+c
b+c
c
b
rb +
r c + r 1 + ra + r 2 ,
= R 1 + ra + r2 =
b+c
b+c
R1 =
R2
donde r1 , r2 representan los errores de redondeo al almacenar los resultados intermedios de las
dos operaciones, mientras que ra , rb , rc son los errores en los datos a, b, c (que contienen tanto
errores inherentes como errores de redondeo). Para la segunda alternativa, se tiene
1 = ra + rb + r1 ,
R
2 = ra + rc + r2 ,
R
ab
ac
3 =
R
R1 +
R2 + r3
ab + ac
ab + ac
b
c
=
(ra + rb + r1 ) +
(ra + rc + r2 ) + r3 .
b+c
b+c
Para obtener las cotas del error, debe conocerse nuestra manera de calcular. Vamos a suponer
que las operaciones se realizan manualmente (base diez), coma flotante con cuatro dgitos para la
parte fraccionaria y redondeando por aproximacion. De este modo, conocemos (por lo comentado
previamente) que la cota del error relativo es r = 21 1014 = 0,0005. Supongamos adem
as que
los datos no contienen error inherente, de manera que ra , rb , rc solo contienen errores de redondeo.
Por tanto, cada una de las cantidades ra , rb , rc , r1 , r2 , r1 , r2 , r3 esta acotada por r. Teniendo esto
en cuenta, podemos acotar
c
b
|rb | +
|R2 |
b + c |rc | + |r1 | + |r2 | + |ra |
b + c
b c
b + c + b + c + 3 r,
b
c
3|
|R
r2 |) + |
r3 |
b + c (|ra | + |rb | + |rc |) + b + c (|ra | + |rc | + |
b c
3
+
+ 1 r.
b + c b + c
15
rn = 13 rn1 ,
n = 1, 2, . . ..
10
3 qn1
qn2 ,
n = 2, 3, . . ..
n = 2, 3, . . ..
17
En la grafica 1.1 se observa una convergencia hacia cero de los errores de la primera sucesi
on.
La figura 1.2 muestra una convergencia de la sucesion de errores correspondientes a la segunda
formula hacia un valor no nulo. Por su parte, la grafica 1.3, asociada a la tercera sucesion, indica
que, a partir de n = 6, los errores empiezan a crecer, siendo inaceptables para valores mayores
de n. Todo ello indica que, para generar la sucesion { 31n }
ormula proporciona el
n=0 , la primera f
mejor procedimiento.
1.5.
Ejercicios
Ejercicio 1. En cada uno de los casos, calcula el error absoluto y el relativo, as como el n
umero
de cifras significativas correctas de la aproximacion.
x = 2,71828182,
x = 98350,
x
= 2,7182.
x
= 98000.
x = 0,000068,
x
= 0,00006.
1/4
x2
e dx
1/4
(1 + x2 +
x4 x6
+ )dx = p.
2!
3!
Determina que tipo de error se presenta aqu y compara el resultado con el valor exacto dado
por el ordenador p = 0,2553074606.
Ejercicio 3. Supongamos que tenemos un ordenador que almacena los n
umeros en base 10 con
tan solo dos dgitos de mantisa. Queremos calcular con esta maquina la menor raz de la ecuaci
on
x2 20x + 1 = 0.
x2 + 1 x, x grande.
cos2 x sin2 x,
x /4.
18
19
Referencias
[1] J.-L.Chabert (Ed.), A History of Algorithms. From the Pebble to the Microchip, Springer,
1999.
[2] P. Henrici, Elements of Numerical Analysis,Wiley, 1964.
[3] A. Huerta, J. Sarrate, A. Rodrguez-Ferran, Metodos numericos. Introduccion, aplicaciones
y programacion, Ediciones UPC, 1999.
[4] J. H. Mathews, K. D. Fink, Metodos numericos con MATLAB, Prentice Hall, 2000.
[5] L. N. Trefethen, Finite Difference and Spectral Methods for Ordinary
and Partial Differential Equations, unpublished text, 1996, available at
http://www.comlab.ox.ac.uk/nick.trefethen/pdetext.html.
[6] L. N. Trefethen, The definition of numerical analysis, SIAM News 25, 6 Nov. 1992; reprinted in Bull. Inst. Maths. and Applics., 1993 and again as an appendix to L. N. Trefethen
and D. A. Bau, III, Numerical Linear Algebra, SIAM, 2000.
[7] J. M. Sanz Serna, Diez Lecciones de Calculo Numerico, Universidad de Valladolid, 1998.
20
Tema 2
2.1.
Introducci
on
ciones para un sistema de 10 ecuaciones). Esta
fue la razon por la que surgieron otros metodos
que reducan el n
umero de operaciones, como el de Gauss (1777-1855), y tambien metodos que
proporcionaban soluciones por aproximacion.
Otra cuestion planteada fue que los coeficientes de estas ecuaciones procedan de medidas,
de mayor o menor precision, y por tanto para obtener la mayor fiabilidad posible, las medidas
(ecuaciones) deberan ser mas numerosas que las incognitas. Era necesario pues encontrar una
forma de sustituir las ecuaciones originales por otro sistema con el mismo n
umero de ecuaciones
que de incognitas y que dieran las aproximaciones a los valores requeridos lo mejor posible. De
entre todos los metodos propuestos, el mas utilizado fue el de mnimos cuadrados de Legendre
(1752-1833) y Gauss. De hecho, Gauss lo utilizo para determinar la orbita de un planeta y
para hacer mas sencillo la obtencion de la solucion de las ecuaciones normales desarroll
o el
metodo de eliminacion que lleva su nombre. Finalmente, otra tarea fue encontrar metodos de
resolucion mas practicos que las reglas de Cramer y la eliminacion gaussiana cuando se necesitaba
manipular sistemas grandes de ecuaciones, en especial los obtenidos usando el metodo de mnimos
cuadrados. Esto dio lugar a la aparicion de los metodos iterativos, desarrollados en el siglo XIX,
que permitan obtener las soluciones con un grado de precision dado. De nuevo Gauss dedujo un
metodo indirecto para resolver un problema de Astronoma (su principal preocupacion). Jacobi
(1804-1851) desarrollo su metodo al considerar sistemas de ecuaciones lineales obtenidas de su
estudio de sistema fsicos sujetos a peque
nas oscilaciones. Por u
ltimo, Seidel (1821-1896), un
alumno de Jacobi, propuso su tecnica iterativa a partir del estudio de un problema de Geodesia,
por lo que parece que no conoca el metodo descrito por Gauss, que era mas popular entre los
astronomos.
Los conocimientos adquiridos sobre la resolucion de sistemas lineales en cursos previos permiten enfocar esta leccion desde un punto de vista numerico y compararla con la siguiente, sobre
procedimientos iterativos. Para algunos conceptos y resultados utilizados (y no explicados) aqu se
recomienda revisar los contenidos del primer curso de algebra lineal.
La estructura de la leccion es como sigue: se presenta primero el concepto de norma matricial
y la idea de acondicionamiento de un sistema. Despues, a una revision general sobre sistemas
lineales y eliminacion gaussiana le seguira la descripcion del algoritmo general, incluyendo la
implementacion del pivotaje. La interpretacion matricial del algoritmo en terminos de operaciones
elementales (analizadas en un apendice final) es la base para la explicacion de la factorizaci
on
LU de una matriz y su aplicacion a la resolucion de un sistema. Finalmente, la aparici
on de
sistemas sobredeterminados permite introducir los problemas de mnimos cuadrados y tratar
su resolucion numerica, a traves, fundamentalmente, de la factorizacion QR. Discutiremos su
aplicacion a problemas de ajuste.
A lo largo de la leccion consideraremos sistemas lineales de ecuaciones de la forma
Ax = b,
A Mmn (C), x Cn , b Cm ,
(2.1)
22
2.2.
2.2.1.
Recordemos que una norma en un espacio vectorial E es una funcion |||| : E R+ verificando
las propiedades
(1) ||v|| 0,
v E con ||v|| = 0 v = 0.
C, v E.
v, w, E.
Figura 2.1: Conjunto de vectores v con ||v||2 = 1 (crculo), ||v|| = 1 (cuadrado), ||v||1 = 1
(rombo).
Las normas para matrices mas utilizadas van ligadas a una norma vectorial de la siguiente
manera. Fijada una norma para vectores || ||, su norma matricial asociada queda definida como
||A|| = sup{||Av||, ||v|| = 1} = sup
v6=0
23
||Av||
.
||v||
(2.2)
Se puede comprobar [4] que (2.2) define una norma en el espacio vectorial de las matrices complejas N N , MN (C), es decir, verifica las propiedades (1)-(3) con matrices en lugar de vectores.
El cociente de (2.2) mide la dilatacion de v debida a la aplicacion de la matriz A. La expresi
on
(2.2) verifica ademas las propiedades
||IN || = 1,
||Av|| ||A||||v||,
||AB|| ||A||||B||,
A MN (C), v CN .
A, B MN (C).
Para las normas matriciales asociadas a las normas vectoriales destacadas anteriormente se pueden
obtener formulas mas explicitas que la definicion general. Sea A = (ajk )N
j,k=1 .
Norma del m
aximo
||A|| = max
1jN
N
X
|ajk |.
k=1
(Se suman los valores absolutos de los elementos de cada fila y se retiene la suma mayor).
En efecto, sea v 6= 0. Por un lado
N
N
X
X
||Av|| = max
ajk vk max
|ajk ||vk |
1jN
1jN
k=1
k=1
max
1jN
N
X
!
|ajk | ||v|| .
k=1
Luego
N
X
||Av||
max
|ajk |,
1jN
||v||
k=1
1jN
N
X
|ajk |.
k=1
Por otro lado, sea J el ndice de la fila en la que se alcanza la mayor suma de los valores absolutos
de los elementos. Se considera el vector y = (y1 , . . . , yN )T con
aJk /|aJk |
si aJk 6= 0
yk =
0
si aJk = 0
donde denota conjugacion. Se tiene entonces que ||y|| = 1, a menos que aJk = 0 para todo
k {1, . . . , N }, en cuyo caso debe ser A = 0, dado que J es la fila con la suma mayor. En estas
condiciones, el resultado es trivial. Para el resto de los casos (aJk 6= 0 para alg
un k {1, . . . , N }),
se tiene
N
N
X
X
||Ay|| = max
ajk yk
aJk yk
1jN
k=1
k=1
!
N
N
N
X
X
X
=
|aJk | = max
|ajk | = max
|ajk | ||y|| .
k=1
1jN
k=1
24
1jN
k=1
Entonces
N
||A||
X
||Ay||
max
|ajk |,
1jN
||y||
k=1
||A||1 = max
1kN
|ajk |.
k=1
(Se suman los valores absolutos de los elementos de cada columna y se retiene la suma mayor).
Norma dos. Para la norma eucldea, es necesario definir lo que se llama radio espectral de una
matriz A MN (C), que es
(A) = max |j |,
j
j R,
||vj ||2 = 1,
1 j N.
N
X
j vj ,
j = hv, vj i, 1 j N,
||v||22
j=1
N
X
|j |2 ,
j=1
se tiene
N
X
j vj A A
j=1
j k vj A Avk =
j,k
N
X
j vj
j=1
j k vj k vk
j,k
2
|j | j (max |j |)
j
|j |2 = (A A)||v||22 .
p
De manera que ||A||2 (A A). Para la otra desigualdad, consideremos un autovalor de
modulo maximo || = (A A) y sea v 6= 0 un autovector asociado. Como
||Av||22 = v A Av = ||v||22 ,
25
||Av||22
= = || = (A A).
||v||22
El radio espectral define la norma dos, pero tambien tiene relacion con el resto de las normas
definidas en la forma (2.2). Si |||| es una norma matricial deducida de una vectorial, A MN (C)
y es un autovalor de A, tomemos un autovector v asociado (v 6= 0 con Av = v). Entonces
||A||
||Av||
||||v||
=
= ||,
||v||
||v||
2.2.2.
4
= ,
3
1+ 8
||A||1 = ,
3
||A||2 = 2.
Acondicionamiento de un sistema
Se dice que una matriz cuadrada invertible A esta mal acondicionada si existe una matriz
B de manera que cambios peque
nos en los elementos de A o B generan cambios grandes en
1
X = A B. Se dice que el sistema Ax = B esta mal acondicionado si A esta mal acondicionada.
En la siguiente leccion comentaremos algo mas sobre el acondicionamiento. Cuando un sistema
esta mal acondicionado, los metodos numericos para su resolucion suelen perder eficacia. Se suele
producir mal acondicionamiento cuando la matriz del sistema es casi singular. Un ejemplo cl
asico
de matrices mal acondicionadas viene dado por las llamadas matrices de Hilbert. La matriz de
Hilbert de orden n es la matriz cuadrada H Mnn (R) de elementos hij = 1/(i + j 1), i, j =
1, . . . , n. Por ejemplo
H(5) =
1/3 1/4 1/5 1/6 1/7 .
1/4 1/5 1/6 1/7 1/8
1/5 1/6 1/7 1/8 1/9
Ejemplo 2.2 [5]. Para el termino independiente b = (1, 0, 0, 0, 0)T y utilizando aritmetica
exacta, la solucion del sistema H(4)x = b es x = (25, 300, 1050, 1400, 630)T . Operando con
siete cifras significativas, la matriz resultante puede considerarse una perturbacion de H(5),
1,0000
0,5000
0,3333
0,2500
0,2000
0,5000000 0,3333333 0,2500000 0,2000000 0,1666667
H(5)
0,3333333 0,2500000 0,2000000 0,1666667 0,1428571 ,
0,2500000 0,2000000 0,1666667 0,1428571 0,1250000
0,2000000 0,1666667 0,1428571 0,1250000 0,1111111
26
(2.3)
Esto significa que el error relativo en x no es mayor que (A) por el error relativo en b, donde
(A) es el llamado n
umero de condicion de A:
(A) = ||A1 ||||A||.
Observemos que el n
umero de condicion depende de la norma matricial elegida. Ademas, si (A)
es peque
no, entonces peque
nos errores relativos en b llevan a peque
nas perturbaciones en x. De
27
este modo, un n
umero de condicion alto significa un mal acondicionamiento del sistema. Hay que
notar que (A) 1 en cualquier norma y para cualquier matriz.
Por ejemplo, la matriz de Hilbert de orden tres
1 1/2 1/3
A = 1/2 1/3 1/4 ,
1/3 1/4 1/5
tiene un n
umero de condicion (A) 524.
1
1/2
A=
1/3
1/4
se tiene que (A) es del orden de 1,6 104 . Las matrices de Hilbert son ejemplos clasicos de
matrices mal acondicionadas.
Al resolver numericamente un sistema de ecuaciones como (2.1), no obtenemos en general
la solucion exacta x sino una aproximacion x
. Nos podemos hacer una idea de la bondad de
tal aproximacion estimando el residuo r = b A
x. Denotemos por e = x x
al error absoluto
cometido. Entonces
Teorema 2.2.1
1 ||r||
||e||
||r||
(A)
.
(A) ||b||
||x||
||b||
Demostraci
on. Notemos primero que Ae = r y que A
x = b, con b = b r. Entonces, la
segunda desigualdad del teorema se deduce de (2.3). Por otro lado,
||r||||x|| = ||Ae||||A1 b|| ||A1 ||||A||||e||||b||,
de donde se deduce la segunda desigualdad.
2.2.3.
Refinamiento iterativo
Cuando se sospecha que la matriz del sistema esta mal acondicionada, existe un procedimiento
que mejora la aproximacion a la solucion del sistema lineal, que se llama refinamiento iterativo.
Consiste en realizar iteraciones en el sistema cuyo lado derecho es el vector residuo para las
aproximaciones sucesivas hasta que se obtiene una precision satisfactoria. Mas concretamente,
sea x(0) una solucion aproximada de (2.1) y r(0) = b Ax(0) . El vector de error e(0) = x x(0)
satisface
Ae(0) = Ax Ax(0) = b Ax(0) = r(0) .
(2.4)
De manera que conocido x(0) , puede calcularse r(0) y posteriormente e(0) resolviendo el sistema
(2.4). Es de esperar que el vector x(1) = x(0) + e(0) sea una aproximacion mas precisa a la soluci
on
de (2.1) que x(0) , pues
x(1) = x(0) + e(0) x e(0) + e(0) .
28
El proceso desarrollado con x(0) puede repetirse para x(1) , generando as aproximaciones mejores
x(2) , x(3) , . . ., a la solucion. El exito de este refinamiento iterativo depende del calculo de los
residuos r(i) con una buena precision, para evitar la cancelacion en la resta.
Ejemplo 2.3. En [4] se considera
420 210
210 140
140 105
105 84
140 105
x1
875
105 84
x2 = 539 ,
84 70 x3 399
70 60
x4
319
2.3.
M
etodo de eliminaci
on gaussiana
Recordamos primero varios aspectos del metodo de eliminacion gaussiana, que necesitaremos
para la descripcion y explicacion del algoritmo asociado. Si bien la eliminacion gaussiana es un
procedimiento programable de discusion y resolucion de un sistema lineal general como (2.1),
esta seccion considerara el caso m = n con A invertible, de modo que (2.1) tiene solucion u
nica.
Escribimos (2.1) como un sistema de ecuaciones
a11 x1 + a12 x2 + + a1n xn
b1
b2
an1 x1 + an2 x2 + + ann xn
(2.5)
bn ,
Recordemos que el metodo de eliminacion gaussiana para este caso consta de dos etapas:
(i) Eliminacion de incognitas y sustitucion progresiva: se pasa del sistema (2.5) a traves de
sistemas equivalentes hasta llegar a un sistema triangular superior.
(ii) Sustitucion regresiva: resolucion escalonada del sistema equivalente de tipo triangular superior, despejando las incognitas desde la u
ltima ecuacion hasta la primera.
29
La primera etapa puede describirse de dos formas, una algebraica y otra matricial. La primera
nos interesara para la descripcion del algoritmo y la segunda al tratar la llamada factorizacion LU
de la matriz A. Desde el punto de vista algebraico, la etapa incluye un proceso de eliminaci
on de
incognitas, con posible intercambio de filas (pivotaje) y una transformacion del termino independiente. Este proceso tiene una estructura algortmica que presentaremos a continuacion. Desde el
punto de vista matricial, el procedimiento consta de la aplicacion sucesiva de operaciones elementales (vease el apendice al final de la leccion) sobre las filas de la matriz ampliada, hasta llevar
el sistema a uno equivalente de tipo triangular superior. Las operaciones elementales constituyen
la version matricial de las operaciones realizadas en el algoritmo, incluyendo el pivotaje.
Ejemplo 2.4. Recordamos lo descrito anteriormente con un ejemplo. Primero aparecen los
sucesivos sistemas equivalentes y despues las correspondientes matrices ampliadas. Los sistemas
equivalentes son:
x1 + 2x2 x3 + x4 = 1
2x1 2x2 + x3 x4 = 1
4x1 10x2 + 5x3 x4 = 1
2x1 14x2 + 4x3 x4 = 0
x1 + 2x2 x3 + x4 = 1
6x2 + 3x3 3x4 = 1
18x2 + 9x3 5x4 = 5
18x2 + 6x3 3x4 = 2
x1 + 2x2 x3 + x4 = 1
6x2 + 3x3 3x4 = 1
4x4 = 2
3x3 + 6x4 = 1
x1 + 2x2 x3 + x4 = 1
6x2 + 3x3 3x4 = 1
3x3 + 6x4 = 1
4x4 = 2
Y la transformacion al sistema triangular superior, desde el
forma siguiente:
1
2
1
2
1 1 | 1
2 2
0 6
1
1
|
1
[A|b] =
4 10 5 1 | 1 0 18
2 14 4 1 | 0
0 18
30
1
3
5
3
| 1
| 1
| 5
| 2
1
0
0
0
2 1
6 3
0
0
0 3
1
3
4
6
|
|
|
|
1
1
0
1
0
2
1
0
2 1
6 3
0 3
0
0
1
3
6
4
|
|
|
|
1
1
.
1
2
2.3.1.
Algoritmo de eliminaci
on gaussiana
31
(k)
(k)
(k)
(k)
(k)
(k)
(k)
(k)
(k)
(k)
(k)
..
.
(k)
(k)
(k)
..
.
(k)
(k)
(k)
(k)
..
.
(k)
(k)
e intercambiamos virtualmente la posicion de las ecuaciones p(l) y p(k), simplemente intercambiando los correspondientes valores en el vector p, quedando este actualizado para el siguiente
pivotaje.
La segunda etapa elimina la incognita xk de las restantes ecuaciones a partir de la p(k)-esima
(la cual, de realizarse el intercambio, estara en la posicion k-esima). Para la eliminacion de xk
en la ecuacion con ndice p(i), i = k + 1, . . . , n, se resta a la ecuacion p(i)-esima lp(i)k veces la
p(k)-esima, donde
(k)
lp(i)k =
ap(i)k
(k)
(2.6)
ap(k)k
lp(i)k =
ap(i)k
(k)
ap(k)k
Fijo i, para j = k + 1, . . . , n
(k+1)
(k+1)
ap(i)j
bp(i)
(k)
(k)
(k)
(k)
Los restantes elementos de las filas p(1), . . . , p(k1) no cambian, mientras que para i = k+1, . . . , n
(k+1)
es ap(i)k = 0 (eliminacion de xk ). Esto completa las otras dos etapas del paso k 7 k + 1.
32
Sustituci
on regresiva
El sistema escalonado final A(n) x = b(n) suele denotarse por U x = c y tendra la forma
up(1)1 x1 + up(1)2 x2 + + up(1)n xn = cp(1)
up(2)2 x2 + + up(2)n xn = cp(1)
..
..
.
.
up(n)n xn = cp(n) .
De manera que la solucion del problema puede obtenerse despejando las incognitas, a partir de
la u
ltima ecuacion hasta la primera.
xn =
xi =
cp(n)
up(n)n
n
X
1
cp(i)
up(i)j xj ,
up(i)i
i = n 1, . . . , 1.
j=i+1
Aspectos de implementaci
on y coste operativo
Ya hemos tratado el intercambio virtual de filas que realiza el ordenador para el pivotaje, a
33
n
X
n
X
n
X
1+
n1
X
n
X
multiplicaciones
k=1 i=k+1
divisiones,
k=1 i=k+1
es decir
n1
X
(n k)2 + (n k) =
k=1
34
n3 n
3
(n k) =
k=1
n2 n
2
j2 =
j=1
N (N + 1)(2N + 1)
.
6
1=n
i=1
n1
X
n
X
divisiones
1=
n1
X
i=1 j=i+1
(n i) =
i=1
n2 n
2
multiplicaciones
es decir, n 2+n operaciones. Observemos entonces que el mayor coste operativo procede de la
transformacion de la matriz A. Cuando esta tiene cierta estructura (veanse los ejercicios al final
de la leccion) puede modificarse el algoritmo para reducir el n
umero de operaciones. A pesar de
3
ello, el coste de la eliminacion (del orden de n /3, siendo n el tama
no del sistema) es muy inferior
al de su metodo predecesor de las reglas de Cramer (recuerdese lo escrito en la introducci
on),
si bien para matrices muy grandes este metodo directo suele sustituirse por metodos iterativos
(tema 4), sobre todo si la matriz tiene una estructura que beneficia el proceso iterativo.
Nota: Al calcular el coste operativo en los algoritmos, suelen aparecer sumas enlazadas con
dos o tres ndices. Si estamos interesados en el termino dominante, es decir, el de mayor potencia
cuando las dimensiones son grandes, este termino puede obtenerse sustituyendo las sumas por
integrales y ajustando los lmites de integracion [7]. Por ejemplo
n X
i
X
n
X
1=
i=1 j=1
i=
i=1
n(n + 1)
n2
2
2
(n
grande),
nZ i
Z
1djdi =
idi =
0
n2
.
2
Hay que observar tambien que el algoritmo puede implementarse sin pivotaje, eliminando
la primera etapa de intercambio de ecuaciones y sustituyendo p(j) por j en todas las formulas
anteriores. Sin embargo, el pivotaje es necesario en bastantes situaciones. Por una parte, no
siempre puede completarse el metodo de eliminacion gaussiana sin pivotaje, dado que es posible
la aparicion de un pivote nulo en el proceso [6]. Pero tambien, desde el punto de vista numerico, un
pivote no nulo pero muy peque
no en relacion con las otras cantidades proporciona multiplicadores
muy grandes que pueden amplificar los errores de la aritmetica finita.
35
Ejemplo 2.5. El siguiente ejemplo se encuentra en [5]. Los valores x1 = x2 = 1 dan la soluci
on
del sistema
1,133x1 + 5,281x2 = 6,414
24,14x1 1,210x2 = 22,93.
Haciendo las operaciones con cuatro cifras significativas, resolvemos el sistema con eliminaci
on
gaussiana sin pivotaje. El u
nico multiplicador a calcular es l21 = 24,14/1,133 = 21,31. Con cuatro
dgitos de precision, el sistema triangular resultante es
1,133x1 + 5,281x2 = 6,414
113,7x2 = 113,8,
el cual, por sustitucion regresiva, da la solucion
x2 =
113,8
= 1,001,
113,7
x1 =
6,414 5,281(1,001)
6,414 5,286
=
= 0,9956.
1,133
1,133
5,338
= 1,000,
5,338
x1 =
22,93 + 1,210(1,000)
= 1,000.
24,14
Por u
ltimo, el pivotaje implementado se denomina parcial, al afectar solo al intercambio de
ecuaciones (o filas de la matriz ampliada del sistema). Puede realizarse tambien lo que se llama
pivotaje total, intercambiando ademas posiciones de las incognitas (o columnas de la matriz
ampliada) [4].
2.3.2.
Interpretaci
on matricial del algoritmo. Factorizaci
on LU
Al principio de la leccion hemos comentado la interpretacion matricial del metodo de eliminacion gaussiana, a traves de las operaciones elementales. Tal interpretacion lleva asociada una
factorizacion de la matriz del sistema, que es la que comentamos en esta seccion. Es conveniente
seguir la explicacion utilizando el apendice 1 del final del tema.
Recordemos que las operaciones elementales sobre las filas de una matriz son:
Intercambiar dos filas.
36
1 0
0
0
0 1
0
0
. .
..
.
.. ..
.
..
1
0
M (k) = 0 0
,
0 0 l
k+1,k
. .
..
..
.. ..
.
.
0 0 ln,k
1
que difiere de la matriz identidad solo en la columna k, la cual tiene la forma
(0, 0, , 0, 1, lk+1,k , , ln,k )T .
Entonces
A(k+1) = M (k) A(k) ,
1
0
0
0
.
..
l21
1
0
0
..
..
..
..
..
..
.
.
.
.
.
.
L= .
.
..
..
..
..
..
..
.
.
.
.
.
l
l
l
1
0
n1,1
n1,2
ln,1
ln,2
n1,n2
ln,n1
Entonces A = LU y b = Lc.
Demostraci
on. Seg
un hemos comentado previamente
U = A(n) = M (n1) A(n1) = = M (n1) M (1) A
c = b(n) = M (n1) b(n1) = = M (n1) M (1) b.
37
1 0
0
0
0 1
0
0
. .
..
..
.. ..
.
1
0
(M (k) )1 = 0 0
,
0 0 l
k+1,k 0
. .
..
.
.. ..
.
..
0 0 ln,k
1
y de nuevo la proposicion 2.6.1 en el apendice 1 nos da el resultado.
La factorizacion LU permite reinterpretar matricialmente las etapas de la eliminaci
on en
gaussiana, en el sentido siguiente:
La eliminacion de incognitas corresponde a llevar a cabo la factorizacion A = LU , con L la
matriz de multiplicadores y U la matriz del sistema triangular superior final. Los elementos
de L y de U pueden almacenarse, como ya mencionamos en la seccion anterior, en la propia
matriz A, de manera que esta, al final, tendra la forma
u11
u12
u1,n1
u1n
l21
u22
u2,n1
u2n
l31
l32
u33
u3,n1
u3n
.
A .
..
..
..
..
..
.
.
.
ln,n1
un,n
El pseudocodigo es entonces (comparese con la seccion anterior)
Algoritmo 2.2.(Continuaci
on). Sustituci
on regresiva
bn bn /ann
Desde i = n 1 hasta 1
desde j = i + 1 hasta n
bi bi bj aij
fin
bi bi /aii
fin
La existencia de la factorizacion LU depende de la posibilidad del desarrollo completo de
la eliminacion gaussiana sin pivotaje. Se puede comprobar que para una matriz cuadrada e
invertible A, existe la factorizacion LU si y solo si todas las submatrices menores principales de
A son invertibles [6].
A la vista de la relacion matricial observada entre los sistemas k y k + 1 de la eliminaci
on
sin pivotaje, puede deducirse la influencia de la incorporacion de este. Hay que tener presente
que eliminar con pivotaje parcial en un sistema es igual que eliminar sin pivotaje en un sistema
obtenido permutando las ecuaciones del original. En este caso aparecen las matrices de permutacion, que son productos de las matrices elementales asociadas a la operacion de permutar
filas.
Ejemplo 2.6. Para el sistema
2x1 3x2 + 14x3 = 4
2x1 + 2x2 3x3 = 5
4x1 + 2x2 2x3 = 10,
39
la eliminacion gaussiana con pivotaje parcial requiere, en su primer paso, intercambiar las
filas 1 y 3, de manera que la primera transformacion del vector de punteros es p = (3, 2, 1). La
inmediata eliminacion de x1 lleva al sistema
4x1 + 2x2 2x3 = 10
x2 x3 = 0
2x2 + 13x3 = 9.
El nuevo sistema requiere un nuevo intercambio entre las filas p(2) y p(3), de modo que ahora
p = (3, 1, 2). Al eliminar x2 , tenemos
4x1 + 2x2 2x3 = 10
2x2 + 13x3 = 9
(2.7)
4,5x3 = 4,5.
Observemos que si denotamos por P la matriz de permutacion que se obtiene intercambiando las
filas de la matriz identidad seg
un p, es decir,
0 0 1
P = 1 0 0,
0 1 0
y si multiplicamos por P a los elementos del sistema original, nos queda
4x1 + 2x2 2x3 = 10
2x1 3x2 + 14x3 = 4
2x1 + 2x2 3x3 = 5,
el cual, resuelto sin pivotaje, nos lleva al sistema equivalente (2.7).
El ejemplo ilustra un resultado general. Si P denota la matriz de permutacion obtenida a
partir de la identidad intercambiando las filas seg
un el u
ltimo valor del puntero de la eliminaci
on
con pivotaje, entonces eliminar con pivotaje parcial el sistema Ax = b equivale a eliminar sin
pivotaje el sistema P Ax = P b. De este modo, los pasos de la eliminacion para este caso son:
Factorizacion P A = LU .
Sustitucion progresiva Lc = P b.
Sustitucion regresiva U x = c.
Pueden identificarse las tres etapas en el programa expuesto en la seccion anterior.
2.4.
2.4.1.
Dos puntos de vista ofrece el planteamiento de esta seccion. Vamos primero a ilustrar la idea
de aproximacion y a centrarnos en la resolucion teorica de los problemas de mnimos cuadrados.
40
Por otro lado, discutiremos tecnicas numericas de resolucion del problema lineal de mnimos
cuadrados para un sistema lineal sobredeterminado, con su aplicacion a los problemas de ajuste.
Los problemas de aproximacion son muy habituales en la Fsica y la Ingeniera. Es destacable
tambien su diversidad, tanto en su planteamiento como en las tecnicas numericas aplicadas para
su resolucion. Sin embargo, buena parte de los problemas de aproximacion pueden establecerse
en un lenguaje matematico com
un. El planteamiento formal sera el siguiente: dado un espacio
real X dotado de una norma, un subconjunto S de X y un elemento x de X, se dice que s S
es una mejor aproximacion a x por elementos de S si
||x s || = nf ||x s||.
sS
La interpretacion del termino mejor aproximacion depende entonces de la norma con la que estemos trabajando, as como de las propiedades del espacio X y del subconjunto S donde buscamos
la mejor aproximacion. Ejemplos tpicos de este planteamiento general son los siguientes.
Aproximaci
on de funciones Un problema importante en Analisis Numerico es la computaci
on
de funciones especiales como las trigonometricas o las hiperbolicas. Ello puede realizarse aproximando la funcion por otra mas sencilla de evaluar. Como ilustracion de esto, podemos buscar
aproximar la funcion sin x en el intervalo [0, /2] por una lnea recta p(x) = ax + b.
Una primera idea consiste en utilizar el desarrollo en serie de Taylor. Podemos escribir
sin x = x
x3 x5
+
.
3!
5!
41
Truncando este desarrollo en el termino lineal, obtenemos nuestra primera aproximacion a sin x:
sin x p1 (x) = x.
La aproximacion esta representada por la lnea de trazo continuo de la figura. Puesto que el
desarrollo de Taylor se toma en torno al cero, la aproximacion es muy buena para valores peque
nos
de x. Sin embargo, para x = /2, el error es aproximadamente 0,57, que la mayora de la gente
considerara inaceptablemente grande.
Otra posibilidad es elegir la recta que pasa por los puntos extremos del grafico de sin x (lnea
discontinua en el dibujo). Esta aproximacion esta dada por
2x
.
(notese que 0,105 0,211/2) obtenemos una aproximacion cuyo error maximo ocurre en tres
puntos: 0, 0,881, /2. Esta aproximacion se muestra en la u
ltima lnea de la figura y es la mejor
posible. Hemos encontrado la solucion al siguiente problema: elegir a y b tal que
p2 (x) =
es mnimo. Notese que este problema de mnimos no se puede resolver por las tecnicas usuales
para el calculo de extremos del calculo infinitesimal, pues F no es una funcion derivable de a, b.
Sistemas sobredeterminados o incompatibles Sea A una matriz real con m filas y n columnas. Tratemos de resolver el sistema A~x = ~b donde ~b es un vector conocido de m componentes
y ~x, desconocido, tiene n componentes. Sea S un subespacio de Rm generado por los n vectores
columna de A. Supongamos que
O bien que n < m, es decir, hay menos incognitas que ecuaciones.
O bien n m, pero entre las n columnas de A no hay m independientes, es decir, hay un
n
umero de incognitas mayor o igual que el de ecuaciones, pero la matriz A no tiene rango
maximo.
Si ~b no esta en S, el sistema es incompatible. Puesto que, cuando ~x recorre Rn , A~x nunca coincide
con ~b, podemos preguntarnos para que ~x es kA~x ~bk lo menor posible, siendo k k una norma
en Rm previamente elegida. Tales ~x juegan el papel de soluciones en un sentido generalizado.
Observemos que su determinacion comprende dos etapas:
(i) Hallar la (o las) mejores aproximaciones ~b a ~b por elementos de la imagen S.
(ii) Resolver el (o los) sistemas A~x = ~b en sentido convencional.
De manera que estos sistemas incorporan, en la etapa (i), un problema de aproximacion.
42
Problemas de ajuste Una aplicacion tpica del apartado anterior viene dada por muchos
problemas de ajuste. Podemos ilustrarlos con el siguiente ejemplo. Supongamos que, cuando se
utiliza una pila de V voltios para cargar un condensador, la diferencia de potencial v entre los
bornes del condensador despues de t segundos satisface una relacion de la forma
v = v(t) = V (V v0 )et/ ,
(2.8)
donde v0 es la diferencia de potencial inicial (en t = 0) y es una constante del circuito. Para
una pila de V = 10 voltios se han obtenido las siguientes mediciones:
t (seg)
0.5
v (volt)
6.36
6.84
7.26
8.22
8.66
8.99
9.43
9.63
A partir de estos datos, nuestro problema es tratar de ajustar los valores de v0 , , as como
utilizar la formula anterior para estimar v cuando t = 6.
No hay ning
un motivo para elegir un dato experimental de la tabla frente a otro, de manera
que ha de contarse con todos ellos y tratar de obtener valores de v0 , y de modo que si
w = (6,36, 6,84, 7,26, 8,22, 8,66, 8,99, 9,43, 9,63)T ,
v = (v(t0 ), v(t1 ), v(t2 ), v(t3 ), v(t4 ), v(t5 ), v(t6 ), v(t7 ))T ,
entonces ||w v|| sea lo menor posible. Si sustituimos directamente los datos en la formula (2.8),
obtenemos un sistema de ocho ecuaciones para dos incognitas, pero no lineal. Una manera de
superar este obstaculo (vease el apendice 2 al final de la leccion) es linealizar el problema. Si
tomamos logaritmos en (2.8), tenemos
t
ln(V v(t)) = ln(V v0 ) ,
(2.9)
ln(10 6,36)
1 0,5
ln(10 6,84)
1 1
ln(10 7,26)
1 2
ln(10 8,22)
1 3
A=
, b = ln(10 8,66) .
1 4
ln(10 8,99)
1 5
ln(10 9,43)
1 7
1
ln(10 9,63)
Una vez resuelto este, se despejan las incognitas = 1/x2 , v0 = V ex1 = 10 ex1 . La resoluci
on
proporciona los parametros = 3,5925 y v0 = 5,7679, as como una estimacion del valor de v en
t = 6 de 9,2034. La figura 2.3 muestra el ajuste proporcionado por la formula (2.8) a los datos
de la tabla.
El ajuste consiste entonces en acercarse a todos los datos a la vez lo mejor posible (en el
sentido dado por la norma) en lugar de pasar por algunos de ellos o por todos. En la practica, el
n
umero de datos a ajustar suele ser mucho mayor y en la mayora de las ocasiones la expresi
on
para v(t) es tambien una hipotesis impuesta a partir de la disposicion de los datos experimentales.
43
Figura 2.3: Datos de la tabla (con estrellas ) y funcion (2.8) dada por el ajuste.
2.4.2.
Los problemas de aproximacion de los que nos ocuparemos se pueden formular en general de la
siguiente forma: dado un espacio vectorial E con un producto interno, ~y E y S = h~v1 , ~v2 , . . . , ~vk i
un subespacio de E, se dice que p~ S es una mejor aproximacion a ~y por elementos de S si
para cada p~ S, es ||~y p~ || ||~y p~||, donde || || es la norma procedente del producto interno
definido en E.
Esto equivale a lo siguiente: p~ = 1~v1 + 2~v2 + + k~vk es una mejor aproximacion a ~y por
elementos de S si 1 , . . . , k son tales que
||~y 1~v1 2~v2 k~vk ||
es mnimo. Tales problemas suelen llamarse de mnimos cuadrados. As, partiendo del planteamiento general, hemos dotado a los elementos con mas propiedades. Ahora el espacio E es vectorial y
eucldeo, la norma procede de un producto interno y el subconjunto de proyeccion S es un subespacio. Veremos despues que estas propiedades son fundamentales para la resolucion de estos
problemas. En el caso particular de ser E = RN y de utilizar la norma eucldea, el problema de mnimos cuadrados tiene la siguiente forma: dados vectores linealmente independientes
~v1 , ~v2 , . . . , ~vk en Rn y un vector ~y en Rn , hay que determinar constantes 1 , . . . , k tales que
||~y 1~v1 2~v2 k~vk ||2
es mnimo. Aqu,
||~u||22 = ~uT ~u =
n
X
i=1
44
u2i ,
h~u, ~v i = ~uT ~v =
ui vi .
i=1
Hay analogos continuos del problema de mnimos cuadrados discreto. Por ejemplo, imaginemos
que queremos aproximar una funcion y(t) por un polinomio de grado k 1 en [1, 1]. Si definimos
~vj (t) = tj1 , j = 1, 2, . . . , k y consideramos la norma
||u(t)||22
u(t)2 dt,
=
1
Z
hu, vi =
u(t)v(t)dt,
1
entonces nuestro problema se reduce otra vez a minimizar ||~y 1~v1 2~v2 k~vk ||2 . En
este caso, se suele tomar como E el conjunto de funciones y : [1, 1] R tales que
Z
y(t)2 dt < ,
n
X
i2 ui vi ,
i=1
||~u||2 =
i2 u2i .
i=1
Entonces, los datos menos precisos perderan influencia en la minimizacion de ||~y 1~v1 2~v2
k~vk ||.
Los problemas continuos tambien pueden ponderarse. Por ejemplo, si deseamos dar m
as importancia a los puntos extremos del intervalo [1, 1], podemos utilizar la norma
2
||u(t)|| =
1
u(t)2
dt.
1 t2
Puesto que la funcion peso 1/ 1 t2 tiende hacia infinito en t = 1, forzara un mejor ajuste
cerca de 1 y 1 que cerca de cero.
45
2.4.3.
Ortogonalidad
Puesto que a partir de un producto interno se puede definir una norma, todo espacio con producto interno es de manera natural un espacio normado, y tiene sentido hablar en el de distancias.
Por otra parte, en los espacios con producto interno hay una serie de conceptos que carecen de
sentido en espacios normados generales. El mas destacado es el de ortogonalidad. Dos vectores
de E se dicen ortogonales si su producto interno es nulo. El concepto de ortogonalidad generaliza
la nocion elemental de perpendicularidad de R2 y R3 y es fundamental para la aproximaci
on en
espacios con producto interno.
]J
J
6
~r
~y JJ
J
J
J
J
J
*
J
(
(( (
(
(
(
(((
~v
(
(
Para ver por que, consideremos el problema de aproximar un vector ~y R3 por una combinacion lineal 1~v1 + 2~v2 de dos vectores linealmente independientes ~v1 , ~v2 . Nuestro criterio para
la mejor aproximacion requiere minimizar la norma del residuo ~r = ~y 1~v1 2~v2 .
La situacion geometrica se ilustra en la figura. Cuando 1 y 2 varan, el vector ~v = 1~v1 +2~v2
se mueve en un subespacio bidimensional representado en la figura por el plano. Las lneas con
trazo mas fino muestran una tpica aproximacion ~v y el correspondiente residuo ~r. La intuici
on
geometrica nos dice que el tama
no de ~r debera hacerse mnimo cuando ~r esta en angulo recto al
plano, como ilustran las lneas de trazo grueso.
En el caso de que la norma derive de un producto interno, puede darse un criterio que nos
permita reconocer si un elemento dado es la mejor aproximacion o no. Siguiendo la intuici
on
geometrica de la figura, tenemos
Sea E un espacio con producto interno, S un subespacio de E, ~y un elemento de E. Si
existe p~ , aproximacion optima a ~y por elementos de S, entonces tal mejor aproximaci
on es
u
nica y satisface
(*)
~y p~ es ortogonal a cada p~ S.
aproximacion.
Es conocido que p~ se llama proyeccion ortogonal de ~y sobre el subespacio S.
Si S es de dimension finita, la relacion (*) da un medio practico para calcular p~ . Tomemos
una base ~v1 , ~v2 , . . . , ~vk de S. Podemos buscar entonces p~ en la forma
p~ =
k
X
i~vi .
i=1
para coeficientes 1 , . . . , k . Para que se verifique (*) se necesita y basta que ~y p~ sea ortogonal
a cada ~vi (al formar estos una base). Por tanto, los i satisfacen
h~y
k
X
i~vi , ~vj i = 0,
j = 1, 2, . . . , k,
j = 1, 2, . . . , k.
i=1
o bien
k
X
(2.10)
i=1
Las ecuaciones (2.10) se llaman ecuaciones normales del problema. La matriz del sistema (h~vi , ~vj i)i,j ,
que es la misma para todas las ~y E, se llama matriz de Gram de los vectores ~vi . El sistema
(2.10) es compatible y determinado, dado que la aproximacion optima existe y es u
nica.
Un caso particular se tiene cuando la base ~vi es ortogonal. Entonces (2.10) toma la forma
||~vj ||2 j = h~y , ~vj i,
j = 1, 2, . . . , k,
con lo que
p~ =
k
X
h~y , ~vj i
j=1
||~vj ||2
~vj .
u(t)2 dt.
=
0
47
es decir
Z
a
Z
a
dt + b
0
Z
1
tdt + b
t2 dt
tdt =
0
1
a + b/2 = 1/3,
0
2
t dt =
A Mm,n (R), ~b Rm , ~x Rn ,
sobredeterminado (m > n, muchas mas ecuaciones que incognitas), con las columnas de A independientes (rango (A) = n) e incompatible, es decir, ~b
/ W = col(A). Por ejemplo,
1 1
5
1 2 15
1 3 x1 = 24 ,
(2.11)
x2
1 4
32
1 5
40
es de este tipo.
Si el sistema A~x = ~b no tiene solucion, entonces ~b no puede ser combinacion de las columnas
de A, es decir, ~b
/ W = col(A). Han de buscarse entonces valores de ~x de modo que, puesto
que A~x no puede ser ~b, este lo mas cerca posible de ~b, en el sentido de que la norma eucldea
||A~x ~b|| sea la menor posible. Tales ~x juegan el papel de soluciones en un sentido generalizado.
Su determinacion comprende dos etapas:
(i) Hallar la mejor aproximacion ~b a ~b por elementos de W (proyeccion ortogonal).
(ii) Resolver el sistema A~x = ~b en sentido convencional. La solucion ~xLS de este sistema se
llama solucion en el sentido mnimos cuadrados, pues por la caracterizacion de la proyecci
on
ortogonal, hace mnima las distancias ||A~x ~b|| cuando ~x recorre Rn .
48
En la practica, estas dos etapas se suelen reunir de la siguiente forma. Denotemos por Ai la
iesima columna de A. Puesto que ~b es la proyeccion ortogonal de ~b sobre W = col(A) y
~b = A~xLS , entonces la condicion de proyeccion ortogonal nos lleva a que ~b ~b = ~b A~xLS debe
ser ortogonal a cada columna de A:
h~b A~xLS , Aj i = 0,
j = 1, 2, . . . , n,
o, equivalentemente
hA~xLS , Aj i = hAj , ~bi,
de modo que las ecuaciones normales en este caso nos llevan a
AT A~xLS = AT ~b.
(2.12)
2.4.4.
Resoluci
on num
erica
a11 a11
Desde j = 2 hasta n
aj1 aj1 /ga11
fin
Desde i = 2 hasta n
desde k = 1 hasta i 1
aii aii a2ik
fin
aii aii
desde j = k + 1 hasta n
desde k = 1 hasta i 1
aji aji ajk
fin
aji aji /aii
fin
fin
Ejemplo 2.9. Para la matriz
1
A = 3
2
3
10
10
50
2
10 ,
21
1 0 0
G = 3 1 0.
2 4 1
La resolucion de las ecuaciones normales (2.12) seguira entonces los siguientes pasos:
1. Calcular la parte triangular inferior de B = AT A (la parte triangular superior se obtiene
por simetra).
2. Calcular d~ = AT ~b.
3. Realizar la factorizacion de Cholesky B = GGT .
~
4. Resolver por sustitucion progresiva G~y = d.
5. Resolver por sustitucion regresiva GT ~xLS = ~y .
El procedimiento requiere del orden de (m + n6 )n2 operaciones.
Factorizaci
on QR de una matriz
El segundo procedimiento se basa en el siguiente resultado.
Cualquier matriz A con columnas independientes puede factorizarse en un producto A =
QR. La matriz Q tiene el mismo tama
no de A y sus columnas son ortonormales (ortogonales
entre s y con norma uno). La matriz R es cuadrada, triangular superior e invertible.
Antes de establecer un ejemplo de la factorizacion QR y su interpretacion algortmica, veamos lo
que significa. Puesto que R es no singular, podemos escribir AS = Q, donde S = R1 . En otras
palabras
snn
o, equivalentemente,
~q1 = s11~v1
~q2 = s12~v1 + s22~v2
l = 1, . . . , k 1.
(b) El nuevo sistema genera el mismo espacio que el de partida. En particular, toda base de E
se puede ortogonalizar.
Demostraci
on. Dados {~v1 , ~v2 , . . . , ~vn } linealmente independientes, el primer vector no cambia:
~e1 = ~v1 .
para obtener ~e2 , escribimos ~e2 = ~v2 12~e1 e imponemos la ortogonalidad entre ~e1 y ~e2 ,
h~e2 , ~e1 i = 0 h~v2 , ~e1 i 12 h~e1 , ~e1 i = 0 12 =
h~v2 , ~e1 i
.
||~e1 ||2
Queda as determinado
~e2 = ~v2
h~v2 , ~e1 i
~e1 ,
||~e1 ||2
que, por construccion, es ortogonal a ~e1 . Ahora, para obtener ~e3 , se escribe ~e3 = ~v3 13~e1 23~e2
y se impone la ortogonalidad con ~e1 y ~e2 ,
h~v3 , ~e1 i
,
||~e1 ||2
h~v3 , ~e2 i
=
.
||~e2 ||2
h~v3 , ~e1 i
h~v3 , ~e2 i
~e1
~e1 .
2
||~e1 ||
||~e2 ||2
El procedimiento permite un paso general. Supongamos que hemos obtenido vectores ~e1 , . . . , ~ek1
ortogonales y queremos obtener ~ek . Entonces escribimos ~ek = ~vk 1k~e1 2k~e2 k1,k~ek1
e imponemos las condiciones de ortogonalidad
h~ek , ~e1 i = 0 1k =
..
.
..
.
h~vk , ~e1 i
,
||~e1 ||2
..
.
h~vk , ~ek1 i
,
||~ek1 ||2
con lo que ~ek queda determinado. Continuamos el proceso hasta obtener tantos vectores como al
principio. .
Fijemonos en que la demostracion proporciona las igualdades
~vk = 1k~e1 + + k1,k~ek1 + ~ek
= 1k ||~e1 ||~q1 + + k1,k ||~ek1 ||~qk1 + ||~ek ||~qk ,
(2.13)
h~v2 , ~e1 i
1
~e1 = ~v2 ~e1 ,
h~e1 , ~e1 i
2
y, operando, es ~e2 = (1/2, 1/2, 1)T . El tercer vector viene de las condiciones
~e3 = ~v3 1,3~e1 2,3~e2
h~e3 , ~e1 i = h~e3 , ~e2 i = 0,
por lo que
~e3 = ~v3
h~v3 , ~e1 i
h~v3 , ~e2 i
1
1
~e1
~e2 = ~v3 ~e1 ~e2 ,
h~e1 , ~e1 i
h~e2 , ~e2 i
2
3
y por tanto ~e3 = (2/3, 2/3, 2/3)T . Los vectores ortonormales finales son entonces
r 1
r 1/2
r 2/3
~e1
1
~e2
2
~e3
3
~q1 =
=
1 , ~q2 =
=
1/2 , ~q3 =
=
2/3 .
||~e1 ||
2
||~e2 ||
3
||~e3 ||
4
0
1
2/3
La factorizacion QR de la matriz A cuyas columnas son ~v1 , ~v2 , ~v3
p
2 p1/2
~v1 ~v2 ~v3 = ~q1 ~q2 ~q3 0
3/2
0
0
53
es
p
1/2
.
p 16
4/3
Algoritmo de Gram-Schmidt
La prueba de la existencia de la factorizacion QR es por tanto constructiva y el algoritmo resultante se llama algoritmo de Gram-Schmidt. Hay que tener en cuenta que, en la implementaci
on,
se calculan directamente los vectores ~qi y los elementos de la matriz R. Por tanto, el paso intermedio de calcular los vectores ortogonales ~ei no se lleva a cabo. De la expresion matricial de la
factorizacion se tiene la igualdad vectorial
~v1 = r1,1 ~q1
~vk =
k1
X
2 k n.
j=1
P
qj ||.
La ortonormalidad de los vectores ~qj nos lleva a que rj,k = h~vk , ~qj i y rk,k = ||~vk k1
j=1 rj,k ~
Una eficiente implementacion exige almacenar los vectores ~qi en los propios ~vi . El pseudoc
odigo
tendra la siguiente forma:
~v1
~vk
1
~vk (~vk contiene el valor de ~qk )
rk,k
Algoritmo modificado
Puesto que el algoritmo de Gram-Schmidt da lugar a la factorizacion QR puede, en principio,
utilizarse para resolver problemas de mnimos cuadrados, como luego veremos. Sin embargo,
este proceso es numericamente inestable: peque
nos errores de redondeo pueden ocasionar que
los vectores calculados no sean ortogonales. Una ligera modificacion del algoritmo evita este
problema. En el algoritmo clasico de Gram-Schmidt, ~q1 se toma m
ultiplo de ~v1 . Entonces, se resta
el m
ultiplo apropiado de ~q1 a ~v2 para obtener un vector que es ortogonal a ~q1 . El procedimiento
modificado de Gram-Schmidt calcula ~q1 como antes. A continuacion, resta m
ultiplos de ~q1 no
(1)
(1)
solo a ~v2 , sino tambien a ~v3 , . . . , ~vn , por lo que los vectores resultantes ~v2 , . . . , ~vn son todos
ortogonales a ~q1 . As, el primer paso del algoritmo modificado de Gram-Schmidt es como sigue
0 = ||~
r1,1
v1 ||
54
~q1 =
v1
0 ~
r1,1
Desde j = 2 hasta n
0 = h~
r1,j
vj , ~q1 i
(1)
~vj
0 ~
= ~vj r1,j
q1
(1)
(1)
(1)
Se puede comprobar facilmente que ~v2 , . . . , ~vn son ortogonales a ~q1 . Puesto que ~v2 es
ortogonal a ~q1 , podemos normalizarlo para obtener ~q2 . A continuacion se resta un m
ultiplo
(1)
(1)
(2)
(2)
apropiado de ~q2 a ~v3 , . . . , ~vn para obtener vectores ~v3 , . . . , ~vn que son ortogonales a ~q2 . Puesto
(2)
(2)
que cada ~v3 , . . . , ~vn es una combinacion lineal de vectores ortogonales a ~q1 , ellos tambien son
ortogonales a ~q1 . Tras k 1 pasos del proceso de Gram-Schmidt modificado tenemos vectores
(k1)
(k1)
. El paso
y vectores ~vk
, . . . , ~vn
que son ortogonales a ~q1 , . . . , ~qk1
ortonormales ~q1 , . . . , ~qk1
(k1)
(k1)
(k1)
(k)
para generar vectores ~vk+1 , . . . , ~vn que son ortogonales a ~q1 , . . . , ~qk . Explcitamente
(k1)
0
= ||~vk
rk,k
||
1 (k1)
vk
0 ~
rk,k
Desde j = k + 1 hasta n
(k1)
0 = h~
vj
, ~qk i
rk,j
~qk =
(k)
~vj
(k1)
= ~vj
0 ~
qk
rk,j
(k)
(k)
, pues cada
Se puede observar que ~vk+1 , . . . , ~vn son ortogonales a ~qk y tambien a ~q1 , . . . , ~qk1
(k)
. Desarrollando estas f
~vj es una combinacion lineal de vectores ortogonales a ~q1 , . . . , ~qk1
ormulas
para k = 1 hasta n, generamos una sucesion de vectores ortonormales ~q1 , . . . , ~qn identicos a los
0
= rj,k . El pseudocodigo es
producidos por el algoritmo de Gram-Schmidt clasico. Tambien rj,k
como sigue:
rn,n = ||~vn ||
~vn
1
rn,n
~vn
55
Se puede comprobar que el coste operativo de ambos algoritmos es del orden de 2n2 m 32 n3
multiplicaciones.
Para la resolucion de,l sistema (2.12), hay que tener presente que si calculamos la factorizaci
on
T
QR de A, entonces, utilizando que las columnas de Q son ortonormales (y, por tanto, Q Q = I)
se tiene
AT A = RT QT QR = RT R, AT ~b = RT QT ~b,
y (2.12) quedara en la forma RT R~x = RT QT ~b. Puesto que R es invertible, el sistema queda
reducido a R~x = QT ~b. La resolucion de (2.12) utilizando factorizacion QR dara entonces los
pasos siguientes:
1. Calcular la factorizacion A = QR.
2. Calcular d~ = QT ~b.
~
3. Resolver por sustitucion regresiva el sistema R~xLS = d.
El mayor coste del procedimiento procede de la factorizacion del primer paso. De este modo,
necesitaramos
Aplicar el algoritmo de Gram-Schmidt para la factorizacion
A = QR 2n2 m
2n3
3
productos.
operaciones.
El coste es el doble que el del metodo de ecuaciones normales si m es mucho mayor que n y aproximadamente el mismo si n = m. Aunque el uso de las ecuaciones normales sea conceptualmente
sencillo y computacionalmente competitivo, puede no ser recomendable por otros aspectos. Uno
de ellos se refiere al hecho de que el n
umero de condicion de AT A puede ser mucho mayor que el
de A. En [4] se considera el ejemplo
1 1 1
1 + 2
1
1
0 0
T
1
A=
1 + 2
1 .
0 0, A A =
1
1
1 + 2
0 0
En A, el parametro 6= 0 evita que su rango sea uno. En AT A, ese parametro es 2 . Cuando
es peque
no, la no singularidad de AT A es computacionalmente mas difcil de identificar que
la de A. Si b = (1, 0, 0, 0)T , se tiene AT b = (1, 1, 1)T y la resolucion de (2.12) para este caso es
1
1
2
2
x1 = x2 = x3 = 3+
no, AT A es computacionalmente singular y
2 = 3 + O( ). Si es muy peque
la resolucion numerica de las ecuaciones normales no funciona. Por contra, para la factorizaci
on
QR, si 2 es despreciable frente a , se tiene que, salvo factores O(2 ), computacionalmente
1 / 2
/ 6
1
1
1
1/ 2 1/ 6
Q=
0 1/ 2 1/ 6 , R = 0 2 /
2 .
0
0
6/2
0
0
2/ 6
2.5.
Ejercicios
Ejercicio 1. Comprueba que si A es una matriz cuadrada definida positiva, entonces la eliminacion gaussiana sin pivotaje aplicada a la matriz A no genera ning
un pivot nulo.
Ejercicio 2. Sistemas tridiagonales. El objetivo es resolver, del modo mas eficiente posible y
sin pivotaje, sistemas lineales con matrices de la forma
d1 s1 0
0
0
0
l1 d2 s2
0
0
0
0 l2 d3
s3
0
0
.. . .
..
..
..
..
..
.
.
.
.
.
.
A= .
0 0 lN 4 dN 3 sN 3
0
0
0 0
0
lN 3 dN 2 sN 2
0
0 0
0
0
lN 2 dN 1 sN 1
0 0
0
0
0
lN 1
dN
Para un sistema A~x = ~b = (b1 , b2 , . . . , bN )T , con A la matriz de arriba, elabora un pseudoc
odigo que para cada una de las tres etapas de la eliminacion en este caso: factorizacion, sustituci
on
progresiva y sustitucion regresiva. Analiza el coste operativo y las necesidades de almacenamiento
del algoritmo. Compara con el algoritmo general.
Ejercicio 3. Resuelve por eliminacion gaussiana (sin pivotaje) el siguiente sistema
2x1 + 2x2 + 2x3 + 2x4 = 12
4x1 + 6x2 + 6x3 + 6x4 = 34
6x1 + 14x2 + 16x3 + 16x4 = 82
2x1 + 2x2 + 12x3 + 14x4 = 42
indicando los valores de los multiplicadores y de los pivotes. Determina la factorizacion LU de la
matriz del sistema.
Ejercicio 4. (i) Sea A una matriz n n conocida e invertible y sea B una matriz n m. Explica
como hallar A1 B sin calcular previamente A1 .
(ii) Se va a calcular y = A1 (Bz + u) + x, donde A y B son matrices n n conocidas y u, x, z
son vectores n-dimensionales tambien conocidos. Explica como disponer los calculos para evitar
invertir A. Compara el costo operativo del proceso sugerido con el obvio de hallar Bz, Bz +
u, A1 (Bz + u) e y.
(iii) Sea A una matriz n n. Especifica un algoritmo que calcule A1 por eliminacion gaussiana
utilizando solo n3 operaciones.
Ejercicio 5. Resuelve por eliminacion gaussiana los siguientes sistemas
4x1 + 5x2 + 6x3 = 28
2x1 7x3 = 29
5x1 8x2 = 64
3x2 13x3 = 50
2x1 6x2 + x3 = 44
4x1 + 8x3 = 4.
57
(2.14)
(2.15)
pi+1 =
ci
,
ai pi + bi
q1 = x0 ,
qi+1 =
di ai qi
.
ai pi + bi
(2.16)
Las formulas (2.15) y (2.16) constituyen el llamado algoritmo de Thomas para sistemas tridiagonales.
Ejercicio 7. Sea A una matriz cuadrada de orden n simetrica y definida positiva. Elabora un
algoritmo que factorice A en la forma A = LDLT donde L es triangular inferior con unos en la
diagonal y D es diagonal con elementos diagonales positivos.
Ejercicio 8. Factorizaci
on de Choleski. Construye un algoritmo para calcular la factorizaci
on
de Cholesky de una matriz cuadrada de orden n simetrica y definida positiva y estudia su coste
operativo.
Ejercicio 9. (i) Demuestra que si aplicamos eliminacion gaussiana sin pivotaje a una matriz
simetrica A, entonces li1 = a1i /a11 .
(ii) A partir del apartado (i), comprueba que si se suprimen la primera fila y la primera
columna de A(2) , se obtiene una matriz simetrica.
(iii) Determina un algoritmo de resolucion de un sistema lineal Ax = b, con A simetrica, con
el menor coste computacional posible.
Ejercicio 10. Resuelve el sistema
2x + 6y = 8,
2x + 6,00001y = 8,00001.
(2.17)
Supongamos que los coeficientes y el segundo miembro se conocen solo aproximadamente (por
ejemplo, a traves de observaciones experimentales) de modo que en un momento dado en lugar
de (2.17) obtenemos el sistema
2x + 6y = 8,
2x + 5,99999y = 8,00002.
(2.18)
1
A = 2
0
2
1
0
0
0,
4
ex f (x)g(x)dx,
utiliza el metodo de Gram-Schmidt para calcular una base ortogonal del espacio de los polinomios
de grado menor o igual que dos con coeficientes reales. Calcula la mejor aproximacion a la funci
on
f (x) = x3 por polinomios ortogonales de segundo grado.
Ejercicio 15. Calcula la factorizacion QR de las matrices
1 1 1
1
1 0 0
,
0
A=
A
=
1 0 2
1
1 1 1
2
1.
4
1 0
A = 0 1,
1 1
y resuelve, en el sentido de mnimos cuadrados, el sistema A~x = ~b = (1, 1, 0)T .
(b) Calcula la factorizacion QR de la matriz
1
0
A=
1
1
1 1
0 1
,
0 0
1 0
59
0
1.0
2
1.6
4
1.4
6
0.6
8
0.2
10
0.8
Ajusta H(t) a las medidas realizadas, empleando el metodo de mnimos cuadrados. Por que no
puede utilizarse directamente la expresion equivalente
2(t t0 )
H(t) = h0 + A sin
?
12
2.6.
Ap
endice 1: Operaciones elementales
Se llama operacion elemental sobre las filas de una matriz a aquella transformacion de la
misma que consiste en llevar a cabo una de las tres manipulaciones siguientes:
(1) Permutar dos filas entre s.
(2) Multiplicar todos los elementos de una fila por un mismo n
umero no nulo.
(3) Sumar a una fila otra fila.
Se sobreentiende que las filas no afectadas por los ndices que definen la operacion elemental
permanecen inalteradas. Fijemonos en que todos los pasos de la eliminacion gaussiana que hemos
dado en los ejemplos anteriores involucran alguno o varios de los tres tipos de operaciones elementales, con ecuaciones en lugar de filas.
En terminos matriciales, la accion de una operacion elemental sobre las filas de una matriz
A equivale a multiplicar por la izquierda la matriz A por una matriz E apropiada asociada a la
operacion. Si A Mm,n (C), la matriz E Mm,m (C) y tiene un aspecto distinto dependiendo
de la operacion elemental. Denotemos por Im a la matriz identidad de m filas y m columnas y,
dados 1 r, s n, por Irs la matriz de Mm,m (C) con todos sus elementos nulos, salvo el que
esta en la posicion (r, s), que vale uno. Entonces:
(1) Si la operacion consiste en permutar las filas r y s (r 6= s), entonces E se obtiene de Im permutando las filas r y s de esta. De este modo, EA es una matriz que se obtiene de A permutando
sus filas r y s. Por ejemplo
3 2
0 0 1
6 7
A = 4 5 , E = 0 1 0 EA = 4 5 .
6 7
1 0 0
3 2
Matricialmente, E = Im + Irs Irr Iss + Isr .
60
2
1
5 8
0 1
,E =
1
0
0
4
EA =
2
4
5 8
0 4
.
2
1
A=
4
5
2 1
1
1
,E = 0
0
1
1
0
1
0
0
2
0 EA = 4
1
0
1
0.
0
1
5
0
A= 4
2
1
5
1
1
1
0 ,E = 0
1
0
0
1
0
1
0
0
EA = 4
1
2
0
5
1
0
0.
1
Toda matriz E Mn,n (K) asociada a una operacion elemental posee inversa. En efecto,
Si se trata de permutar dos filas, E es su propia inversa.
Si se trata de multiplicar una fila dada por 6= 0, entonces la inversa de E es la matriz
asociada con la multiplicacion de la misma fila por 1/.
Si se trata de sumar a la fila r la fila s (r 6= s), la inversa es F EF , donde F es la matriz
asociada con la operacion elemental: multiplicacion de la fila s por (1).
Por ejemplo:
0 0 1
E = 0 1 0 E 1 = E,
1 0 0
1 0 0
1 0 0
E = 0 2 0 E 1 = 0 1/2 0 ,
0 0 1
0 0 1
1 0 1
1 0 0
1
1
E= 0 1 0
E = 0 1 0
0
0 0 1
0 0 1
0
0
1
0
1
1
0
0
1
0
0
1
0
0
0 .
1
Por otra parte, para el tipo de matrices asociadas a la eliminacion gaussiana, es decir, de la
1 = I I .
forma Ers = I +Irs , 6= 0 (sumar a la fila r la fila s multiplicada por ), entonces Ers
rs
61
En efecto, si w = Ers v, entonces todas las componentes de w coinciden con las de v salvo la resima, que vale wr = vr + vs . Esto implica que vj = wj , j 6= r y vr = wr ws , es decir
1 w = (I I )w.
v = Ers
rs
Otra propiedad relevante de estas matrices es la siguiente:
Proposici
on 2.6.1 Si r > s, r0 > s0 , r < r0 y Ers = I + Irs , Er0 s0 = I + Ir0 s0 , entonces
Ers Er0 s0 = I + Irs + Ir0 s0 .
(2.19)
Demostraci
on. Sean w = Er0 s0 v, z = Ers w. Entonces
wj = vj ,
j 6= r0 ,
j 6= r,
zr = wr + ws .
zj = wj ,
zr = wr + ws = vr + vs ,
lo que demuestra (2.19).
Estas dos propiedades son utilizadas en la demostracion del teorema 2.3.1 de factorizaci
on
LU .
2.7.
Ap
endice 2. Algunas herramientas estadsticas
En este segundo apendice nos proponemos dar una forma de cuantificar los errores producidos
en la aproximacion por mnimos cuadrados en los problemas de ajuste mas utilizados: la regresi
on
lineal y la polinomial. Vease por ejemplo [1] y referencias incluidas en el.
2.7.1.
Fundamentos estadsticos
10.08
10.11
10.09
10.10
10.12
9.88
10.14
10.02
9.80
10.21
10.19
9.79
9.69
10.05
9.78
10.04
9.98
10.02
9.97
10.04
x
=
1X
xi .
n
i=1
n
X
(xi x
)2 .
(2.20)
i=1
Por lo tanto, si las medidas individuales se dispersan muy lejos de la media, S, y por tanto sx
crecera. Si se agrupan muy cerca de la media, entonces la desviacion estandar sera peque
na. La
dispersion tambien se puede representar por el cuadrado de la desviacion estandar, a la que se
llama varianza
S
s2x =
.
n1
Una medida estadstica final que tiene utilidad en la cuantificacion de la dispersion de los datos
es el coeficiente de variacion c.v.: es el cociente de la desviacion estandar de la media. Como tal,
proporciona una media normalizada de la dispersion. A menudo se multiplica esta cantidad por
100, de tal manera que se puede expresar en forma porcentual
c.v. =
sx
100 %.
x
Estas medidas estan relacionadas con los errores producidos en los experimentos. Volviendo al
ejemplo del principio de la seccion, en la tabla siguiente se calculan la media y la desviaci
on
estandar de la muestra experimental proporcionada por cada persona.
Persona
A
10.12
Media
10.10
Desviacion
0.016
10.08
10.11
10.09
10.10
9.88
10.14
10.02
9.80
10.21
10.01
0.17
10.19
9.79
9.69
10.05
9.78
9.90
0.21
10.04
9.98
10.02
9.97
10.04
10.01
0.033
Los resultados obtenidos por A tienen dos caractersticas importantes. Primero, todos est
an
muy proximos; todos caen entre 10.08 y 10.12. La segunda caracterstica es que todos son demasiado altos (recordemos que la respuesta correcta es 10.00). Han surgido entonces dos tipos
63
de errores distintos en el experimento de A. En primer lugar hay errores aleatorios o accidentales, los cuales provocan que los resultados individuales caigan a ambos lados del valor medio
(10.10). Se dice que estos errores afectan a la precision o reproducibilidad de un experimento. En
el caso de A, queda claro que los errores aleatorios son peque
nos y, por lo tanto, decimos que los
resultados son precisos. La desviacion de los experimentos con respecto a la media es peque
na.
Esto se traduce en un valor peque
no para la desviacion estandar de la muestra. Sin embargo,
tambien existen errores sistematicos, los cuales provocan que todos los resultados sean err
oneos
en el mismo sentido (en este caso todos son demasiado altos). Los errores sistematicos afectan a
la exactitud, es decir, a la proximidad al valor verdadero. Esto se traduce en una media alta.
Analizamos brevemente los resultados obtenidos por las otras personas. B ha obtenido resultados en contraposicion directa a los de A. El promedio de los cinco datos (10.01) esta muy pr
oximo
al valor verdadero, de manera que podemos caracterizar los datos como exactos, sin errores sistematicos sustanciales. Sin embargo, la variedad de los resultados es grande (alta desviaci
on)
lo que indica una precision insatisfactoria y la presencia de errores aleatorios destacables. La
comparacion de estos resultados con los obtenidos por A muestra con claridad que los dos tipos
de errores pueden ocurrir independientemente uno del otro. Esta conclusion se refuerza con los
datos de C y D. El trabajo de C no es ni preciso ni exacto, mientras que D ha logrado a la vez
resultados precisos y exactos.
Distribuci
on de los errores
Aunque la distribucion estandar proporciona una medida de la dispersion de un conjunto de
resultados alrededor del valor medio, no indica la forma en que estan distribuidos los mismos.
Para aclarar esto, se requiere un gran n
umero de mediciones como las expuestas en la siguiente
tabla, en la que se proporcionan los resultados de 50 determinaciones repetidas de concentraci
on
de una cierta sustancia en una muestra concreta de agua, en miligramos por mililitro.
0.51
0.51
0.51
0.50
0.51
0.49
0.52
0.53
0.50
0.47
0.51
0.52
0.53
0.48
0.49
0.50
0.52
0.49
0.49
0.50
0.49
0.48
0.46
0.49
0.49
0.48
0.49
0.49
0.51
0.47
0.51
0.51
0.51
0.48
0.50
0.47
0.50
0.51
0.49
0.48
0.51
0.50
0.50
0.53
0.52
0.52
0.50
0.50
0.51
0.51
Concentracion
0.46
Frecuencia
1
0.47
0.48
0.49
10
0.50
10
0.51
13
0.52
0.53
tiene una desviacion estandar simbolizada por y el valor de la desviacion estandar de la muestra,
sx , nos proporciona una estimacion de .
exp((x )2 /2 2 )
,
2
cuya forma se muestra en la figura 2.5. La curva es simetrica respecto a y cuanto mayor sea el
65
Un analisis un poco mas detallado muestra que, cualesquiera que sean los valores de y ,
aproximadamente el 68 % de los valores de la poblacion caen en el intervalo (, +); cerca del
95 % de los valores se ubican dentro de (2, +2) y casi el 99.7 % de los valores se encuentra en
(3, +3). Esto significa que si las concentraciones de las sustancias que aparecen en la tabla
se distribuyen normalmente, cerca del 68 % caeran en el intervalo (0,483, 0,517); alrededor del
95 % en el intervalo (0,467, 0,533) y el 99.7 % en (0,450, 0,550). De hecho, 33 de los 50 resultados
(66 %) caen entre 0,483 y 0,517; 49 (98 %) entre 0,467 y 0,533, y todos los resultados se ubican
entre 0,450 y 0,550, de manera que la concordancia con la teora es bastante satisfactoria.
Aunque no se pueda demostrar que las mediciones repetidas de cualquier cantidad analtica
siempre van a estar distribuidas normalmente, las pruebas indican que por lo general esta hip
otesis
esta cerca de la verdad.
Intervalos de confianza
Ya hemos visto que la media de una muestra de mediciones nos proporciona una estimaci
on
del valor verdadero , de la magnitud que se quiere medir. Sin embargo, ya que las mediciones
individuales estan distribuidas en torno al valor verdadero con una dispersion que depende de
la precision, es poco probable que la media de la muestra sea exactamente igual al valor verdadero. Por esta razon, es mas u
til proporcionar un intervalo de valores que contenga casi con
seguridad el valor verdadero. La amplitud de este intervalo depende de dos factores. El primero
es la precision de las mediciones individuales, las cuales dependen a su vez de la varianza de
66
la poblacion. El segundo es el n
umero de mediciones de la muestra. El solo hecho de repetir
mediciones implica que se tiene mas confianza en la media de varios valores que en uno solo.
Si consideramos los resultados de la concentracion de sustancia como 10 muestras en las que
cada una de ellas contiene cinco resultados, podemos ver que las mediciones de las muestras de
este tama
no estan distribuidas alrededor de . Si tomamos cada columna como una muestra,
las medias son 0,506, 0,504, 0,502, 0,496, 0,502, 0,492, 0,506, 0,504, 0,500, 0,486. Estas mediciones
se encuentran mas agrupadas entre s que las mediciones originales. De la misma forma que estas
eran muestras de una poblacion infinita de posibles mediciones, estas medias son una muestra
de las posibles medias de muestras de cinco mediciones de la poblacion global. La distribuci
on
de estas mediciones muestrales se denomina distribucion muestral de la media; su media es la
misma que la de la poblacion original y su desviacion estandar se denomina error estandar de
la media (e.e.m.). Existe una relacion matematica entre el e.e.m. y la desviacion estandar, ,
de la distribucion de mediciones individuales, la cual es independiente de la forma en que esten
distribuidas. Si n es el tama
no de la muestra, la relacion es
e.e.m. = / n.
1,96(/ n) < x
< + 1,96(/ n)
Sin embargo, en la practica disponemos de una muestra, de media conocida y buscamos un
intervalo para . Puede reordenarse la ecuacion y expresarse
+ 1,96(/ n).
x
1,96(/ n) < < x
Esta ecuacion proporciona el intervalo de confianza al 95 % de la media. De manera similar, el
intervalo al 97 % esta dado por
x
2,58(/ n) < < x
+ 2,58(/ n).
En el caso de nuestra muestra, tenemos x
= 0,500, n = 50. Lo u
nico que se desconoce es .
Para muestras grandes, sx proporciona una estimacion de y puede sustituirse por esta. As, el
intervalo de confianza al 95 % para la concentracion de sustancia es
67
x
t(/ n) < < x
+ t(/ n).
La variable t representa una nueva distribucion, conocida como distribucion de Student. Depende
tanto de n 1, que se conoce como n
umero de grados de libertad, como del grado de confianza
requerido. Para valores grandes de n, esta distribucion puede aproximarse por la normal. Los
valores de t, al igual que los de la normal, estan tabulados en funcion de los grados de libertad
y el grado de confianza. El termino grados de libertad se refiere al n
umero de desviaciones
independientes (xi x
) que se utilizan al calcular sx . En este caso, dicho n
umero es n 1,
porque
cuando se conocen n 1 desviaciones, la u
ltima se puede deducir si se utiliza el resultado
P
(xi x
) = 0. En la siguiente tabla se recogen algunos valores de t. En ella puede apreciarse que
para muestras de tama
no superior a 50, los valores de t son muy proximos a los de 1.96 y 2.58
utilizados en las ecuaciones anteriores. Esto confirma la validez de la hipotesis utilizada antes al
calcular los lmites de confianza para la concentracion de sustancia.
2.7.2.
Grados de libertad
1
2
3
4
5
10
20
30
50
100
95 %
12.71
4.30
3.18
2.78
2.57
2.23
2.09
2.04
2.01
1.98
99 %
63.66
9.92
5.84
4.60
4.03
3.17
2.85
2.75
2.68
2.63
Regresi
on lineal
x1
y1
x2
y2
xn
yn
n
X
(yi a bxi )2 .
i=1
68
(2.21)
Los parametros a y b se obtienen calculando la solucion mnimos cuadrados del sistema sobredeterminado
y1
1 x1
y2
1 x2 a
.
=
..
... .
..
. b
1 xn
yn
Cuantificaci
on del error en la regresi
on lineal
Resuelto el problema anterior, se puede derivar una serie de propiedades adicionales del ajuste.
Notese la similitud entre (2.20) y (2.21). En el primer caso, los residuos representaban la diferencia
entre los datos y una aproximacion simple de la medida de la tendencia central: la media. En
la ecuacion (2.21), los residuos representan el cuadrado de la distancia vertical entre los datos y
otra medida de la tendencia central: la lnea recta. La analoga se puede extender mas para casos
en donde: 1) la dispersion de los puntos alrededor de la recta son de magnitud similar a lo largo
del rango entero de los datos y 2) la distribucion de estos puntos alrededor de la lnea es normal.
Una desviacion estandar de la lnea de regresion se puede determinar como
r
Sr
sy/x =
.
n2
Esta ecuacion difiere de la de la desviacion estandar sx de una serie de mediciones repetidas en que
las desviaciones son sustituidas por los residuos y en que ahora, en un calculo de regresion lineal,
el n
umero de grados de libertad es n 2, lo que refleja el hecho obvio de que para representar una
recta solo se necesitan dos puntos. La variable sy/x se llama error estandar de la aproximaci
on. Al
igual que la desviacion estandar, el error estandar de la aproximacion cuantifica la dispersi
on de
los datos. Sin embargo, sy/x cuantifica la dispersion alrededor de la lnea de regresion, mientras
que la desviacion estandar sx cuantifica la dispersion alrededor de la media.
Los conceptos anteriores se pueden emplear para medir la eficiencia del ajuste. Esto es particularmente u
til en la comparacion de varias regresiones. Los errores aleatorios en los valores
(2.22)
x
,
b+x
(2.23)
70
Funcion y = f (x)
Linealizacion Y = Ax + B
Cambios
y = A/x + B
y = A/x + B
X = 1/x, Y = y
y = C1 (xy) +
y = D/(x + C)
D
C
X = xy, Y = y
C = A1 , D = B
A
y=
1
Ax+B
1
y
= Ax + B
X = x, Y = 1/y
y=
x
Ax+B
1
y
= A x1 + B
X = 1/x, Y = 1/y
y = A ln(x) + B
y = A ln(x) + B
X = ln(x), Y = y
y = CeAx
ln(y) = Ax + ln(C)
X = x, Y = ln(y)
C = eB
y = CxA
X = ln(x), Y = ln(y)
C = eB
y = (Ax + B)2
y = CxeDx
y=
L
1+CeAx
y 1/2 = Ax + B
ln xy = Dx + ln(C)
ln
L
y
1 = Ax + ln(C)
X = x, Y = y 1/2
X = x, Y = ln xy
C = eB , D = A
X = x, Y = ln Ly 1
C = eB
71
2.7.3.
Regresi
on polinomial
n
X
i=1
Para obtener los coeficientes del polinomio optimo, ya sabemos que se tiene que calcular la
solucion mnimos cuadrados del sistema sobredeterminado
a0
1 x
x21 xk1
y1
1
a
1 x2 x22 xk2 1 y2
.
a2 = . .
..
..
..
..
..
.
.
.
. .. ..
.
yn
1 xn x2n xkn
ak
As como en la regresion lineal, el error en la regresion polinomial se puede cuantificar mediante
el error estandar de la aproximacion
s
Sr
sy/x =
,
n (k + 1)
donde k es el orden del polinomio. Dado que se usaron k +1 coeficientes a0 , . . . , ak , se han perdido
k + 1 grados de libertad. Ademas del error estandar, se puede calcular tambien el coeficiente de
correlacion en la regresion polinomial
r2 =
donde
Sv =
Sv Sr
,
Sv
n
X
(yi y)2 .
i=1
72
Referencias
73
Tema 3
Interpolaci
on polin
omica
1. Introduccion.
2. Interpolacion de Lagrange.
3. Interpolacion de Hermite.
4. Interpolacion trigonometrica.
5. Ejercicios.
6. Apendice. Series de Fourier.
3.1.
Introducci
on
Este tema aborda el problema de la interpolacion polinomica. En muchas situaciones, uno debe
plantearse el estudio de un modelo del que solo conoce datos tabulados, procedentes de mediciones
experimentales. Una manera de estudiar el fenomeno, con vistas a trabajar sobre el (mas all
a de
los datos) es adaptarlo, de alguna forma, a una representacion continua del problema, buscando
funciones de variable continua que tengan alguna relacion. Aqu trataremos de la interpolaci
on
polinomica, que es una forma de representar los datos experimentales. La interpolacion busca un
polinomio que pase por todos ellos. La eleccion de los polinomios como representacion se debe
sobre todo a que estan dentro de las funciones mas sencillas de manejar.
La interpolacion es una idea que ya se utilizaba en la Antig
uedad, sobre todo para obtener tablas matematicas con las que facilitar los calculos. Por ejemplo, en Astronoma, tablas
trigonometricas como las de Ptolomeo en el siglo II d. C., han jugado un importante papel
historico, al igual que la introduccion de tablas logartmicas para facilitar los calculos. Todas ellas utilizaban alg
un tipo de interpolacion entre los datos, para conocer los valores de las funciones
requeridas (las tablas de Ptolomeo, por ejemplo, utilizaban interpolacion lineal, mientras que hay
tratados
arabes de Astronoma, del siglo V d. C, que interpolan senos de angulos con formulas
de orden dos). A medida que las tablas requeran mayor precision, los metodos de interpolaci
on
precisaban de un mayor desarrollo. Las formulas clasicas se establecieron entre el final del siglo
XVII y principios del XIX (Newton, Lagrange, Cauchy).
74
Hora
10
12
14
16
18
Grados
12
18
21
19
15
Una situacion ilustrativa del papel de la interpolacion viene dada en el siguiente ejemplo. La
tabla 3.1 representa la medicion, a ciertas horas del da, de la temperatura ambiente en cierta zona
de cierta ciudad. Buscamos estimar la temperatura a las 13 horas. Una manera de hacerlo consiste
en representar esta como una funci
on del tiempo t, considerado como una variable continua, a
traves de un polinomio que interpole los datos de la tabla. El valor de ese polinomio en el tiempo
t = 13 nos proporcionara una estimacion de la temperatura a la hora solicitada.
Se pueden plantear diferentes problemas de interpolacion, en funcion del tipo de polinomios
que se utilizan. En este tema trataremos tres: el problema de interpolacion de Lagrange, el
problema de interpolacion de Hermite y la interpolacion trigonometrica (quedan fuera problemas
como la interpolacion segmentaria). Esta u
ltima es de especial interes en la teora de se
nales. En
ella se cambia la interpolacion con polinomios en potencias de x (utilizada en los dos primeros
problemas) por los polinomios trigonometricos, en los que las funciones de representacion son
sinusoides de frecuencias m
ultiplos de una conocida.
El tratamiento de los tres problemas sera similar: analizaremos primero la resolucion as como
la representacion del polinomio obtenido; luego trataremos aspectos de implementacion numerica, para terminar discutiendo el nivel de aproximacion del polinomio a la funcion interpolada,
mediante el analisis del error.
3.2.
Interpolaci
on de Lagrange
j = 0, . . . , N.
(3.1)
x0
x1
xN
y0
y1
yN
75
= y0 ,
aN xN
1
= y1 ,
..
.
a0 + a1 x1 + +
..
.
a0 + a1 xN + +
aN xN
N
(3.2)
= yN ,
donde, recordemos, los xj son todos distintos. Esto implica que el sistema (3.2) tiene soluci
on
u
nica y, por tanto, hay un u
nico polinomio PN de grado menor o igual que N y que resuelve el
problema 1.
3.2.1.
Representaci
on del polinomio
Lj (xi ) = 0, i 6= j.
j = 0, . . . , N,
(3.3)
Para completar (3.3), necesitamos expresiones explcitas de los Lj . Observemos que, por definicion, Lj (x) tiene N ceros x0 , x1 , . . . , xj1 , xj+1 , . . . , xN , de modo que ha de ser de la forma
Lj (x) = j (x x0 ) (x xj1 )(x xj+1 ) (x xN ).
76
Y la imposicion de la u
ltima propiedad Lj (xj ) = 1 determina la constante j ,
j =
1
,
(xj x0 ) (xj xj1 )(xj xj+1 ) (xj xN )
quedando entonces
Lj (x) =
Ejemplo 2.1. Polinomio de grado menor o igual que seis que interpola la tabla 3.1, escrito en
forma de Lagrange. El polinomio tiene la forma
P (x) = 7L0 (x) + 9L1 (x) + 12L2 (x) + 18L3 (x)
= +21L4 (x) + 19L5 (x) + 15L6 (x).
Los polinomios de Lagrange para este caso son
L0 (x) =
=
L1 (x) =
=
L2 (x) =
=
L3 (x) =
=
L4 (x) =
=
L5 (x) =
=
L6 (x) =
=
7680
(x 6)(x 8)(x 12)(x 14)(x 16)(x 18)
(10 6)(10 8)(10 12)(10 14)(10 16)(10 18)
(x 6)(x 8)(x 12)(x 14)(x 16)(x 18)
.
3072
(x 6)(x 8)(x 10)(x 14)(x 16)(x 18)
(12 6)(12 8)(12 10)(12 14)(12 16)(12 18)
(x 6)(x 8)(x 10)(x 14)(x 16)(x 18)
.
2304
(x 6)(x 8)(x 10)(x 12)(x 16)(x 18)
(14 6)(14 8)(14 10)(14 12)(14 16)(14 18)
(x 6)(x 8)(x 10)(x 12)(x 16)(x 18)
.
3072
(x 6)(x 8)(x 10)(x 12)(x 14)(x 18)
(16 6)(16 8)(16 10)(16 12)(16 14)(16 18)
(x 6)(x 8)(x 10)(x 12)(x 14)(x 18)
.
7680
(x 6)(x 8)(x 10)(x 12)(x 14)(x 16)
(18 6)(18 8)(18 10)(18 12)(18 14)(18 16)
(x 6)(x 8)(x 10)(x 12)(x 14)(x 16)
.
46080
77
Una ventaja de la representacion (3.3) frente a la dada por el sistema (3.2) es que no necesita
la resolucion de un sistema lineal de tipo Vandermonde para obtener los coeficientes. Tiene, sin
embargo y al menos, un inconveniente. La representacion (3.3) no es apropiada para la evaluaci
on
del polinomio, un problema que s resolva la representacion en potencia de x, a traves del
algoritmo de Horner, presentado en el tema anterior.
Ejemplo 2.2. Se puede contar el n
umero de operaciones necesarias para evaluar PN en la
forma (3.3). El coste mayor procede de la computacion de los polinomios de Lagrange. Para cada
j = 0, . . . , N , la evaluacion de Lj (x) requiere 2(N 1) multiplicaciones y una division, de manera
que en total tenemos
N
X
2(N 1) = 2(N 2 1)
multiplicaciones,
j=0
N
X
1 = N +1
divisiones,
j=0
yN PN 1 (xN )
.
(xN x0 ) (xN xN 1 )
(3.4)
(3.5)
yj Pj1 (xj )
(xj x0 ) (xj xj1 )
yj (c0 + c1 (xj x0 ) + c2 (xj x0 )(xj x1 ) + + cj1 (xj x0 ) (xj xj1 ))
=
,
(xj x0 ) (xj xj1 )
j = 1, . . . N.
(3.6)
La expresion (3.5) se denomina forma de Newton del polinomio interpolador. Los coeficientes
(3.6) se llaman diferencias divididas. Una de las ventajas de la forma (3.5) es que estas diferencias
divididas pueden calcularse a traves de un algoritmo, que pasamos a describir.
Escribamos cj = f [x0 , . . . , xj ], dado que, por (3.6), cj depende solo de los j + 1 primeros
nodos de interpolacion. El siguiente resultado proporciona el algoritmo para el calculo de los
coeficientes.
Teorema 3.2.1 (i) El valor de f [x0 , . . . , xj ] no depende del orden de los nodos, es decir, la
funci
on f [x0 , . . . , xj ] es simetrica con respecto a los nodos x0 , . . . , xj .
(ii) Se tiene adem
as que f [x0 ] = y0 y
f [x0 , . . . , xj ] =
j = 1, . . . , N.
(3.7)
y0
y1
y2
..
.
xN
yN
f [x0 , x1 ]
f [x1 , x2 ]
..
.
f [x0 , x1 , x2 ]
..
.
..
.
f [xN 1 , xN ] f [xN 2 , xN 1 , xN ] f [x0 , . . . , xN ]
Los elementos de la diagonal de la tabla son los coeficientes del polinomio interpolador PN en la
3,60
1,80
1,20
0,90
0,72
3,60
1,80 1,80
1,20 0,60
0,90 0,30
0,72 0,18
0,60+1,80
31
0,30+0,60
42
0,18+0,30
53
= 0,60
= 0,15
= 0,06
0,150,60
41
0,06+0,15
52
= 0,15
= 0,03
0,03+0,15
51
= 0,03
El polinomio interpolador es
P4 (x) = 3,60 1,80(x 1) + 0,60(x 1)(x 2) 0,15(x 1)(x 2)(x 3)
+0,03(x 1)(x 2)(x 3)(x 4).
80
3.2.2.
Si los datos de interpolacion corresponden a valores de una funcion en los nodos del problema
yj = f (xj ),
j = 0, . . . , N,
(x x0 ) (x xN ) N +1)
f
().
(N + 1)!
(3.8)
Demostraci
on. Sea x [a, b]. Si x coincide con alg
un xj , entonces ambos lados de la igualdad
(3.8) se anulan para cualquier . Si x 6= xj , j = 0, . . . , N , se define
M=
f (x) PN (x)
,
(x x0 ) (x xN )
y la funcion
F (t) = f (t) PN (t) M (t x0 ) (t xN ).
Hay que observar que
1. F N +1) (t) = f N +1) (t) (N + 1)!M .
2. F (xj ) = 0,
j = 0, . . . , N,
F (x) = 0.
Por la segunda propiedad, aplicando reiteradamente el teorema de Rolle ([3], por ejemplo), se
obtiene un punto con las condiciones del enunciado tal que F N +1) () = 0. Por la primera
propiedad, se tiene entonces f N +1) () = (N + 1)!M , lo cual, teniendo en cuenta la expresi
on de
M , prueba la igualdad (3.8).
Comentarios
1. Dado que el punto de la formula (3.8) es en general desconocido, esta se utiliza como cota
de error: suponiendo que f N +1) esta acotada en [a, b], para cada x [a, b], se tiene
|f (x) PN (x)|
|x x0 | |x xN |
max f N +1) ()
(N + 1)!
[a,b]
2. Las diferencias divididas proporcionan otra expresion para el error de interpolacion. Dado
x [a, b], sea Q el polinomio de grado menor o igual que N + 1 que interpola a f en
x0 , . . . , xN , x. Entonces
Q(t) = PN (t) + f [x0 , . . . , xN , x](t x0 ) (t xN ).
De modo que, evaluando en x
f (x) PN (x) = f [x0 , . . . , xN , x](x x0 ) (x xN ).
(3.9)
Comparando las expresiones (3.8) y (3.9) se puede observar que, salvo la constante factorial,
las diferencias divididas proporcionan aproximaciones a las derivadas.
81
3.3.
Interpolaci
on de Hermite
r
X
mj ,
j=0
y tal que las siguientes derivadas de f existen, se busca determinar un polinomio PN (x), de grado
menor o igual que N , tal que para j = 0, . . . r
PN (xj ) = f (xj ),
m )
N)
P0 (x) = 2,
P1 (x) = P0 (x) + c1 (x 1), P10 (1) = c1 = 1,
P2 (x) = P1 (x) + c2 (x 1)2 = 2 + (x 1) + c2 (x 1)2 ,
P2 (3) = 2 + 2 + 4c2 = 1 c2 = 3/4.
3
P3 = P2 (x) + c3 (x 1)2 (x 3) = 2 + (x 1) (x 1)2 + c3 (x 1)2 (x 3),
4
P30 (3) = 1 3 + 4c3 = 2 c3 = 1.
Luego
3
P3 (x) = 2 + (x 1) (x 1)2 + (x 1)2 (x 3).
4
Notese que para leer las condiciones de interpolacion, puede tomarse cualquier orden. En lo
u
nico que afecta esto es en la base en la que se expresa el polinomio interpolador. Una ordenaci
on
que utilizaremos mas adelante viene dada al leer las condiciones por columnas, en lugar de por
filas, como se ha hecho en la demostracion (es decir, imponiendo primero las condiciones sobre f ,
luego las de la derivada, luego las de la derivada segunda, etc). Esto significa expresar el polinomio
83
en la base
1, (x x0 ), (x x0 )(x x1 ), . . . , (x x0 ) (x xr ),
(x x0 )2 (x x1 ) (x xr ), . . . , (x x0 )2 (x x1 )2 (x xr )2 ,
3.3.1.
f j) (x0 )
.
j!
(iii) Si entre los argumentos de f [x0 , . . . , xj ] hay dos distintos (que podemos suponer, por (i), el
primero y el u
ltimo) entonces se tiene la formula
f [x0 , . . . , xj ] =
j = 1, . . . , N.
Ejemplo 2.5. Tabla de diferencias divididas y coeficientes del polinomio del ejemplo 2.4. La
tabla se escribe
1 2
1 2
3 1
3 1
1
1/2
2
3/4
5/4
1/2
1/2
2 1/2
3/4
5/4
tenemos los coeficientes del polinomio en la base 1, (x1), (x1)(x3), (x1)2 (x3). (Justifica
la formacion de la cuarta columna de la tabla).
3.3.2.
Estimaci
on del error
3.4.
(x x0 )N +1 N +1)
f
().
(N + 1)!
Interpolaci
on trigonom
etrica
85
3.4.1.
Problema de interpolaci
on trigonom
etrica
Los elementos del nuevo problema de interpolacion son los llamados polinomios trigonometricos. Un polinomio trigonometrico DM de grado M sobre un intervalo [a, b] es una funcion de la
forma
M
X
2im
DM (x) =
cm exp
(x a) ,
(3.10)
ba
m=M
(3.11)
m=1
k = 0, 1, 2 . . . , M,
bk = i(ck ck ),
k = 1, 2 . . . , M.
Hay que observar que DM es una funcion periodica de perodo L = b a. Las expresiones
(3.10)-(3.11) tienen su justificacion en el contexto de la representacion de se
nales. Es razonable
pensar, siguiendo la teora de series de Fourier (vease el apendice de esta leccion) que si (3.10) o
(3.11) va a representar una se
nal, contenga frecuencias tanto positivas como negativas.
Problema 3. Sean L = b a > 0 y N 1 entero. Sobre el intervalo [a, b] se toman los nodos
L
, j = 0, . . . , N 1. Conocidos yj C, j = 0, . . . , N 1, se busca
equiespaciados xj = a + j N
determinar un polinomio trigonometrico de grado N que interpole la tabla 3.3.
x
x0
x1
xN 1
y0
y1
yN 1
La hipotesis de que los nodos xj esten equiespaciados no es restrictiva en la practica. Tpicamente, los datos yj son valores muestreados de una se
nal, cuya expresion en terminos de una
funcion de variable continua es normalmente desconocida, de forma que la b
usqueda de una representacion continua de la se
nal con muestra yj admite la libertad de eleccion de la posici
on de
los nodos que, por sencillez desde el punto de vista matematico y tambien por alg
un que otro
argumento practico, pueden suponerse equiespaciados.
Como en los casos anteriores, el primer punto a discutir es la existencia de solucion u
nica del
problema. Para ello, debemos comprobar el siguiente resultado.
2
N N ,
Lema 3.4.1 Sea = N = exp 2i
= cos 2
N
N i sin N , N 1 y la matriz FN C
N
N
2i(k 1)(l 1)
(k1)(l1)
FN =
= exp
.
(3.12)
N
k,l=1
k,l=1
Entonces:
86
(i) FN es simetrica.
N
(ii) las
columnas de FN forman un sistema ortogonal de C . La norma de cada columna es
N.
1
N FN ,
donde FN = F N = F N .
Demostraci
on. Es claro, de la propia definicion, que FNT = FN , de manera que se tiene (i).
Para (ii), tomemos dos columnas ek , el , 0 k, l N 1 de FN . Teniendo en cuenta que
= 1 ,
entonces
hek , el i =
N
1
X
k j
l j
( ) (
) =
N
1
X
j=0
(kl)j ,
j=0
FNT F N = N IN = F N FN ,
con IN la matriz identidad de tama
no N . Como FN es simetrica, entonces
1
1
FN
F N = IN =
F N FN ,
N
N
lo que implica que FN es invertible y que FN1 =
1
N FN.
N
1
X
ck exp
k=0
2ik
(x a) .
L
(3.13)
Demostraci
on. La comprobacion es constructiva y proporciona la expresion para los coeficientes, que mas tarde utilizaremos. Imponiendo los datos de interpolacion de la tabla 3.3 a una
expresion de la forma (3.13), tenemos
yj = PN 1 (xj ) =
N
1
X
ck exp
k=0
87
2ik
(xj a) .
L
Sea = exp
2i
N
= cos
2
N
i sin
2
N
. Entonces
2ik
2ikj
exp
=
jk ,
(xj a) = exp
L
N
de manera que
yj = PN 1 (xj ) =
N
1
X
ck
jk ,
0 j N 1.
(3.14)
k=0
(3.15)
con la matriz FN dada por (3.12). Para comprobar el teorema, bastara entonces con ver que el
sistema (3.15) tiene solucion u
nica. Por el lema 3.4.1, FN es regular y simetrica, luego FN = FN
es tambien invertible y ademas, tambien por el lema 3.4.1, (FN )1 = N1 FN . Por tanto, el sistema
(3.15) tiene solucion u
nica y esta dada por
c = (FN )1 y =
1
FN y.
N
Ya hemos comentado que es razonable pensar que la se
nal a representar contenga frecuencias
tanto positivas como negativas, de manera que la representacion del polinomio trigonometrico
dada por el teorema 3.4.2 parece incompleta, pues el polinomio carece de frecuencias negativas, y
por tanto da la impresion de que solo aproximara apropiadamente a aquellas funciones cuyos coeficientes de Fourier de ndice negativo sean nulos, es decir, funciones que no contengan frecuencias
negativas. Ello fuerza a buscar representaciones del polinomio interpolador mas convenientes.
Supongamos, por razones que mas adelante justificaremos, que N = 2M es par. Entonces,
una representacion mas adecuada del polinomio es
M
1
X
2ik
TM (x) =
ck exp
(x a) .
(3.16)
L
k=M
1
f f tshif t(FN y),
N
o equivalentemente
f f tshif t(c) =
1
(FN y).
N
88
Demostraci
on. Si f f tshif t(c) = (
c0 , . . . , cN 1 )T entonces
ck
si 0 k M 1
ck =
ckN
si M k N 1
Notemos ademas que para cada par de ndice k, j, se tiene
(
)kj = (
)(N k)j .
Imponemos ahora las condiciones de interpolacion. Para 0 j N 1,
yj
M
1
X
= TM (xj ) =
kj
ck (
)
k=M
M
1
X
k=0
M
1
X
kj
ck (
)
ck (
)kj +
k=0
M
X
M
1
X
kj
ck (
)
1
X
k=0
kj
ck (
)
k=1
N
1
X
clN (
)lj =
l=M
ck (
)kj
k=M
M
1
X
k=0
N
1
X
kj
ck (
)
M
X
ck (
)(N k)j
k=1
ck (
)kj .
k=0
TM
(x)
M
1
X
1
2ik
2iM
ck exp
cM exp
(x a) +
(x a)
2
L
L
k=M +1
1
2iM
cM exp
(x a) ,
2
L
(3.17)
que divide el primer termino de (3.16) en dos aportaciones con frecuencias opuestas. La representacion (3.17) tiene presente que el coeficiente de Fourier de una se
nal real asociado a una
frecuencia es conjugado del asociado a la frecuencia opuesta (vease el apendice). Si los datos
(x) = Re(T (x)) es real, a diferencia de T (x), lo que
de interpolacion son reales, entonces TM
M
M
parece dar una representacion mas fiel de las se
nales reales.
Ejemplo 2.6. Para la funcion f (x) = 1/(5 4 cos x), la figura 3.1 muestra las tres formas del
polinomio trigonometrico interpolador de f en los puntos xj = 2j/8, j = 0, . . . , N 1 = 7. Dado
que la funcion es real y tiene frecuencias positivas y negativas, la tercera forma (3.17) parece la
mas apropiada.
3.4.2.
La matriz FN dada por (3.12), que aparece en el calculo de los coeficientes del polinomio
trigonometrico interpolador, obtenido en la seccion anterior y en cualquiera de sus formas, es de
especial importancia en el analisis numerico y en la teora de la representacion de se
nales. La
matriz define una transformacion lineal FN : CN CN , conocida como transformada de Fourier
89
Figura 3.1: Se
nal muestreada y las tres formas del polinomio trigonometrico interpolador.
N
1
X
ul e
2ik
l
N
l=0
N
1
X
kl
ul N
,
k = 0, . . . , N 1.
l=0
fk =
M
1
X
2jk
u2j N
+
j=0
= (FM u
M
1
X
(2j+1)k
u2j+1 N
j=0
(p)
)k +
k
N
(FM u(i) )k
M
1
X
jk
k
u2j M
+ N
j=0
(p)
fk
90
k (i)
N
fk ,
M
1
X
jk
u2j+1 M
j=0
k = 0, . . . , N 1,
con
u(p) = (u0 , u2 , . . . , uN 2 )T , u(i) = (u1 , u3 , . . . , uN 1 )T .
k+M
k , las f
Por u
ltimo, como N
= N
ormulas anteriores pueden dividirse en dos,
(p)
(i)
(p)
(i)
k
fk = fk + N
fk ,
fk+M
k
= fk N
fk ,
(3.18)
k = 0, . . . , M 1.
(3.19)
De este modo, el calculo de una DFT de N componentes se reduce al calculo de dos DFT de N/2
componentes. Este proceso puede repetirse si, a su vez, M es par. De este modo, si N es una
potencia de dos, digamos N = 2r , entonces el procedimiento puede aplicarse r veces.
El algoritmo consiste entonces en descomponer la se
nal original de N componentes, tras r
pasos, en N se
nales de una componente, calcular la DFT de estas se
nales (algo inmediato, pues
como F1 = 1, la DFT de una se
nal de una componente es ella misma). Finalmente, se agrupan
estas DFT unidimensionales en la DFT de la se
nal original, en r pasos, utilizando las formulas
(3.18) con M = 20 , 21 , . . . , 2r1 . Vamos a analizar cada una de estas tres etapas.
Descomposici
on de la se
nal o bit reverso
La figura 3.2 ilustra, con N = 8 = 23 , el proceso de reordenacion de las componentes de la
se
nal original para su descomposicion en N se
nales de una componente. La se
nal de dieciseis
componentes se descompone en r = 4 etapas. La primera rompe la se
nal en dos se
nales de ocho
componentes, la segunda en cuatro se
nales de cuatro componentes, la tercera en ocho de dos y,
finalmente, la cuarta, en dieciseis se
nales de una componente. La filosofa es la siguiente: cada
vez que una se
nal debe dividirse en dos, se separa en sus componentes pares y sus componentes
impares.
u0
u1
u2
u3
u4
u5
u6
u7
@
R
@
u0
u2
u4
u4
@
R
@
u0
u1
@
R
@
u0
u6
u2
u6
u1
u2
u5
@
R
@
u6
u1
u5
u7
@
R
@
@
R
@
u4
u3
u3
u7
@
R
@
u5
u3
u7
Decimal
0
1
2
3
4
5
6
7
Binaria
000
001
010
011
100
101
110
111
Binario(reverso)
000
100
010
110
001
101
011
111
Decimal
0
4
2
6
1
5
3
7
@
@
@
@
@
@
@
@
@u
a bW
a + bW
92
TM (x ) =
M
1
X
ck z0k = zM cM + cM +1 z0 + cM 1 z02M 1 .
k=M
El algoritmo FFT ha ido admitiendo diversas variantes, vease, por ejemplo, [2].
3.4.3.
Cotas de error
X
k=
ck exp
2ik
(x a) ,
L
93
x [a, b]
(3.20)
ck =
1
L
Z
a
2ik
f (x) exp
(x a) dx,
L
L = b a,
y supongamos que la serie converge absolutamente a f (x) para todo x [a, b]. Supongamos
tambien que N = 2M y consideremos el polinomio trigonometrico (3.13) de coeficientes dados
por c = (1/N )FN y.
Lema 3.4.4 (Aliasing). En las condiciones antes citadas, para cada n = 0, . . . , N 1,
cn =
cn+pN .
(3.21)
p=
Demostraci
on. Observemos primero que la serie de (3.21) tiene sentido al suponer que la
serie de Fourier de f converge absolutamente. Vamos a ver que el polinomio trigonometrico QN 1
de coeficientes
dn =
cn+pN ,
n = 0, . . . , N 1,
p=
X
X
2in
(xk a)
QN 1 (xk ) =
cn+pN exp
L
n=0 p=
!
N
1
X
X
2ink
=
cn+pN exp
N
p=
n=0
!
N
1
X
X
2ik(n + pN )
=
cn+pN exp
N
n=0 p=
X
2ikj
=
cj exp
= f (xk ).
N
j=
Las dos u
ltimas igualdades se deben a que cualquier ndice j Z puede escribirse en la forma
j = n + pN para cierto n = 0, . . . , N 1 y cierto p Z (la pen
ultima) y a la convergencia de la
serie de Fourier hacia f en [a, b] (la u
ltima).
La formula (3.21) se conoce como formula del aliasing (superposicion). Demuestra que todas
las frecuencias n + pN, p Z, estan representadas en la frecuencia n del polinomio trigonometrico
de coeficientes dados por la transformada de Fourier discreta. Si los coeficientes de Fourier cj
mayores corresponden a frecuencias bajas en modulo (algo tpico en muchas se
nales en la practica)
on
entonces cj es una buena aproximacion a cj para j = 0, . . . , N2 1, pero cj es buena aproximaci
a cjN para j = N/2, . . . , N 1. Ello explica la conveniencia de TM (con frecuencias positivas y
negativas) para aproximar mejor a la funcion que PN 1 (con solo frecuencias positivas).
94
Error en la interpolaci
on trigonom
etrica
Dada una funcion g : [a, b] C vamos a considerar las normas
||g|| =
sup |g(x)|
axb
1/2
Z
||g||2 =
|g(x)| dx
||g||2 L||g|| .
(3.22)
Buscamos ahora, como en los anteriores casos de interpolacion, estimar la diferencia |f (x)TM (x)|
para puntos x [a, b] (o bien en R, por la periodicidad). Se dice que la interpolacion es de orden
> 0 si existe C > 0 tal que
|f (x) TM (x)| C/N ,
x R, N = 2M = 2, 4, . . .
x R, N = 2M = 2, 4, . . .
||f TM ||2 2 L
|
ck |
|k|M
Demostraci
on. Escribimos la diferencia f TM = f SM + SM TM , donde SM es el
polinomio trigonometrico truncado obtenido a partir de la serie de Fourier de f ,
M
1
X
2ik
SM (x) =
ck exp
(x a) .
L
k=M
P
Dada la convergencia, el termino f SM puede acotarse por |k|M |
ck |. Para el segundo termino,
se utiliza la formula del aliasing,
|SM (x) TM (x)|
M
1
X
k=M pN,p6=0
|
ck+pN |
|
ck |,
|k|M
X
X
C
dx
|
ck | 2
2C
(3.23)
m+2
m+2
k
M 1 x
k=M
|k|M
xm1
1
|
.
= 2C
=O
m + 1 M 1
N m+1
Esto, junto con el lema 3.4.5, prueba las dos primeras igualdades. Notemos, por otro lado, que los
coeficientes de Fourier de TM son ck para k = M, . . . , M 1 y nulos para los restantes valores
de k. Entonces, utilizando la formula del aliasing y (3.23),
X
X
1
|
ck ck | =
|
ck+pN |
|
cj | = O
.
N m+1
pN,p6=0
|j|M
Cuando la funcion f es analtica, se obtiene una aproximacion de orden exponencial, algo que
no se consigue con ninguna de las interpolaciones anteriores.
96
4C
2 B
L
NB
L
1e
Adem
as, para todo ndice k,
|
ck ck |
3.5.
2C
exp(N B/L).
1 exp(2B/L)
Ejercicios
Ejercicio 1. Comprueba que el valor del determinante de la matriz del sistema (3.1) es
xj ).
i<j (xi
Ejercicio 2. Comprueba que los polinomios de Lagrange L0 , . . . , LN forman una base del espacio
de polinomios de grado menor o igual que N y coeficientes complejos. Comprueba ademas que
se verifican las igualdades
N
X
Li (x) = 1,
i=0
N
X
xki Li (x) = xk , k = 0, . . . , N.
i=0
-1
2
0
1
1
2
2
-7
3
10
A
nade un termino a p de modo que el polinomio que resulte interpole a la tabla entera.
Ejercicio 4. Comprueba que si f es un polinomio de grado menor o igual que N , entonces,
cualquier diferencia dividida de f de orden mayor o igual que N + 1 es nula.
Ejercicio 5. Se desea construir una tabla para la funcion f (x) = ex con abscisas uniformemente
espaciadas en el intervalo [0, 1] de longitud h. Determine cuanto debe valer h si se desea que el
error cometido al hacer interpolacion lineal en dicha tabla sea menor que 5E 05.
Ejercicio 6. Prueba que si se toma cualquier conjunto de 23 nodos en el intervalo [1, 1] y se
interpola all la funcion f (x) = cosh x con un polinomio P (x) de grado a lo sumo 22, el error
relativo que se comete no es superior a 5 1016 .
Ejercicio 7. Un polinomio P (x), cuyo grado se desconoce, admite valores dados por la siguiente
tabla:
97
x
0 1 2
3
P(x) 4 9 15 18
Si todas las diferencias divididas de orden cuatro son iguales a uno, calcula el coeficiente de
x3 en P (x).
Ejercicio 8. Se pretende evaluar la funcion
Z
1
J0 (x) =
cos(x sin t)dt,
0
en el punto 2,15. Se dispone de la siguiente tabla de valores:
2,0
2,1
2,2
x
J0 (x) 0,2239 0,1666 0,1104
Calcula el polinomio interpolador de grado dos de la tabla, emplealo para aproximar el valor
J0 (2,15) y acota el error obtenido.
Ejercicio 9. (i) Escribe el polinomio de Taylor de grado dos de f (x) = ex en x0 = 0.
(ii) Escribe en potencias de x 1 el polinomio obtenido en el apartado (i).
(iii) calcula el polinomio de Taylor de grado dos de f (x) en x0 = 1. Coincide con el anterior?
Ejercicio 10. Sea P (x) un polinomio. Recordemos que se dice que P tiene una raz (un cero) de
multiplicidad m en x = a si P (x) = (x a)m Q(x) con Q(a) 6= 0. Comprueba que P (x) tiene una
raz de multiplicidad m si y solo si P (a) = P 0 (a) = = P m1) (a) = 0.
Ejercicio 11. Calcula el polinomio de Hermite que interpola los datos
x
P (x)
P 0 (x)
0
2
9
1
4
4
2
44
Ejercicio 12. Construye un polinomio de grado menor o igual que tres tal que P (1) = 1, P 0 (1) =
2, P 00 (1) = 4, P (2) = 5 en las bases
{1, (x 1), (x 1)2 , (x 1)3 }.
{1, (x 1), (x 1)(x 2), (x 1)2 (x 2)}.
{1, (x 2), (x 1)(x 2), (x 2)(x 1)2 }.
Ejercicio 13. Construye un algoritmo para calcular los coeficientes del polinomio interpolador de
Hermite y evaluarlo en un punto para el caso mj = 1, j = 0, . . . , r. Estima el coste computacional
en terminos del grado N del polinomio.
Ejercicio 14. Consideremos el polinomio trigonometrico, en el intervalo [a, b],
P (x) =
N
1
X
k=0
ck exp
2ik(x a)
L
,
98
L = b a,
jL
N ,j
= 0, . . . , N 1. Comprueba que
Ejercicio 15. Se consideran 2N + 1 nodos x0 = 0 < x1 < < x2N < 2 y valores complejos
y0 , y1 , . . . , y2N . Comprueba que existe un u
nico polinomio de la forma
N
a0 X
P (x) =
+
(aj cos(jx) + bj sin(jx)) ,
2
j=1
x0
x1
x2 xN 2 xN 1
xN 1
x0
x1 xN 3 xN 2
X = xN 2 xN 1 x0 xN 4 xN 3
x1
x2
x3 xN 1
x0
entonces
~x ~y = X T ~y .
Comprueba las siguientes propiedades de la convolucion:
1. x y=y x.
2. Si f = FN x, g = FN y, h = FN (x y), entonces hj = fj gj , j = 0, . . . N 1.
Ejercicio 20. Sea f : R R una funcion L-periodica, cj el j-esimo coeficiente de Fourier
continuo de f . Supongamos que
|
cj |
1
,
(1 + j 2 )5
j = 0, 1, 2, . . . .
M
1
X
ck exp(2ikx/T ).
k=M
99
(1) f : [0, 2] R,
(2) f : [0, 2] R,
1
.
(1 + cos2 (x))
3
f (x) =
.
5 + 4 cos x
f (x) =
(3.24)
3.6.
Ap
endice. Series de Fourier
Se incluye aqu un apendice sobre la teora de series de Fourier. Esta teora tiene su importancia cuando las funciones a las que se aplica representan se
nales analogicas. Desde el punto de
vista matematico, su aplicacion se extiende a la resolucion de ecuaciones diferenciales. La primera
referencia a la teora se debe a Newton. En el a
no 1672 publico un artculo describiendo el experimento de descomposicion de la luz blanca en los colores del arco iris, utilizando un prisma de
cristal (vease la figura 3.6). Newton empleo el termino espectro para describir las bandas de luz
de colores producidas por el prisma. Luego coloco otro prisma invertido con respecto al primero
y mostro que los colores volvan a mezclarse para producir luz blanca. Insertando una ranura
100
opaca entre ambos prismas y bloqueando la incidencia de uno o mas colores sobre el segundo
de ellos mostro que la luz producida a la salida de este ya no era blanca. Newton dedujo que lo
que estaba produciendo el primer prisma era la descomposicion de la luz en ciertas componentes
elementales, y, el segundo, la recomposicion de estas.
Todos estos experimentos, desde un punto de vista matematico, conllevan la separacion de una
se
nal (en este caso la luz) en componentes elementales de tipo sinusoidal, de distintas frecuencias,
mediante el prisma. As, en el experimento de Newton, el primer prisma realiza el analisis en
frecuencia de la se
nal, obteniendo su espectro. El segundo prisma realiza la sntesis de la se
nal a
partir de su espectro.
El Analisis de Fourier trata entonces de representar funciones mediante una serie de senos y
cosenos, o funciones sinusoidales, llamada serie de Fourier. La finalidad de esta representaci
on se
encuentra al interpretar la funcion representada como una se
nal. El Analisis de Fourier identifica
una serie de elementos en la se
nal que en otros tipos de representaciones no estan tan claros. Estos
elementos son muy u
tiles no solo para identificar la se
nal, sino tambien para almacenarla (algo
muy importante) transmitirla, manipularla, recuperarla, etc. En primer lugar, determinaremos
la forma de calcular la serie de Fourier de una funcion dada. Luego, como para todo tipo de
aproximacion, desde el punto de vista matematico, hay que aclarar en que sentido la serie de
Fourier representa a la funcion. Terminaremos con algunos comentarios sobre las series de Fourier
de senos y de cosenos.
3.6.1.
Definici
on y c
alculo de una serie de Fourier
Aunque las definiciones que vienen a continuacion pueden extenderse a un tipo mas amplio,
nos limitaremos a considerar aqu funciones reales f (x) definidas en un intervalo [a, b] y que
son continuas a trozos. Esto significa que en el intervalo [a, b] existe un n
umero finito de puntos
x1 < . . . < xN 1 tal que si x0 = a, xN = b, f es continua en cada subintervalo (xi1 , xi ), i =
101
Z
f (x)dx = 2
f (x)dx.
0
f (x)dx = 0.
L
Definici
on. dada una funcion f Cp ([a, b]), la serie trigonometrica
a0 X
2
2
+
ak cos
k(x a) + bk sin
k(x a)
,
2
L
L
(3.25)
k=1
donde L = b a y
ak =
bk =
2
f (x) cos
k(x a) dx,
L
a
Z
2 b
2
f (x) sin
k(x a) dx,
L a
L
2
L
k = 0, 1, . . . ,
(3.26)
k = 1, 2, . . . ,
(3.27)
102
(1) f : [, ] R, f (x) = x.
Z
1
a0 =
xdx = 0,
Z
1
x cos(k(x + ))dx = 0,
ak =
Z
1
1 x cos(k(x + )) sin(k(x + ))
2
bk =
x sin(k(x + ))dx =
| = .
+
2
k
k
k
Entonces
sin(2(x + )) sin(3(x + ))
1
+
+ + sin(k(x + )) + ).
2
3
k
1 + x si 1 x 0
(2) f : [1, 1] R, f (x) =
.
1x
si 0 x 1
Z 1
Z 1
a0 =
f (x)dx = 2
(1 x)dx = 1,
f (x) = x 2(sin(x + )
1
1
ak =
1
1
2
(1 (1)k ),
(k)2
bk =
1
De modo que
1
4
f (x) 2
2
1
1
cos((x + 1)) + cos(3(x + 1)) + +
cos((2k + 1)(x + 1)) + .
9
(2k + 1)2
(3) f : [1, 1] R,
f (x) = |x|.
Z 1
Z
a0 =
f (x)dx = 2
1
Z 1
xdx = 1,
ak =
1
1
2
2 k2
((1)k 1),
Z
bk =
Por tanto
1
1
cos((x + 1)) + cos(3(x + 1)) + +
cos((2k + 1)(x + 1)) + .
9
(2k + 1)2
0 si < x < 0
(4) f : [, ] R, f (x) =
.
x
si 0 < x <
Z
Z
1
1
a0 =
f (x)dx =
xdx = ,
0
2
Z
1
1
ak =
f (x) cos(k(x + ))dx = 2 ((1)k 1),
k
Z
1
1
bk =
f (x) sin(k(x + ))dx = .
k
4
1
f (x) + 2
2
103
Figura 3.8: Graficas de las funciones y sumas parciales de la serie de Fourier para los ejemplos
(1) y (2).
De modo que
1
1
f (x)
cos(x + ) + cos(3(x + )) + +
cos((2k + 1)(x + )) +
9
(2k + 1)2
sin(2(x + )) sin(3(x + ))
1
+(sin(x + ) +
+
+ + sin(k(x + )) + ).
2
3
k
2
+
4
ck ei L k(xa) ,
(3.28)
k=
con
1
ck =
L
(3.29)
se llama serie de Fourier compleja de f . Las relaciones entre los coeficientes de la anterior forma
trigonometrica (3.25) y la forma compleja (3.29) se obtienen a partir de las formulas de Euler
cos(t) =
eit + eit
,
2
sin(t) =
eit eit
,
2i
t R.
De manera que
c0 =
a0
,
2
ck =
ak ibk
,
2
ck =
104
ak + ibk
, k = 1, 2, . . . ,
2
(3.30)
Figura 3.9: Graficas de las funciones y sumas parciales de la serie de Fourier para los ejemplos
(3) y (4).
o bien
ak = ck + ck ,
bk = i(ck ck ),
k = 0, 1, 2, . . .
(3.31)
k = 1, 2, . . .
(3.32)
Hay que observar que la representacion (3.25) tiene un sumatorio que recorre los n
umeros naturales (k = 1, 2, . . .) mientras que (3.28) contiene un sumatorio que recorre todos los n
umeros
enteros (k = . . . , 2, 1, 0, 1, 2, . . .). Esto es esencial para que ambas representaciones sean equivalentes.
Ejemplos.
(1) f : [, ] R,
mente, se tiene
c0 = 0,
i
i
ck =
bk =
,
2
k
i
i
ck =
bk = , k = 1, 2, . . .
2
k
Entonces
f (x) = x
X
i
(1)k+1 eik(x+) .
k
k=
105
(2) f : [, ] R,
ck =
si x 0
.
si 0 < x
Z
1/2
1
ik(x+)
0
f (x)e
dx =
(1)k
2
f (x) =
0
1
ik
si
si
si
k=0
k 6= 0 es par
k es impar
Entonces
f (x)
i X e(2k+1)i(x+)
1
+
.
2
2k + 1
k=
Notas:
La existencia de los coeficientes de Fourier de una funcion f (x) no justifica la convergencia
de la serie. Por ello es usual la notacion .
La cuestion acerca de en que sentido la serie de Fourier representa a la funcion se tratar
a en
las siguientes secciones.
Las componentes de una serie de Fourier son funciones periodicas de perodo L. De manera
que, si la serie tiene sentido, la funcion resultante debe ser periodica de perodo L. Por tanto,
no podemos esperar que la serie converja a f (x) para todo x, sino mas bien a la extensi
on
periodica de f . Si estamos interesados en representar por medio de una serie de Fourier
una funcion f definida en [a, b], podemos extender dicha funcion a todo R de manera que
la extension sea una funcion periodica de perodo L. Para ello, basta repetir la grafica de
f en intervalos de longitud L. En la figura 3.8 se visualiza la extension 2 periodica de
f : [, ] R, f (x) = x.
3.6.2.
106
Aproximaci
on mnimos cuadrados
Utilizando la formula de Euler
ei = cos + i sin ,
podemos expresar un polinomio trigonometrico de grado N en la forma
TN (x) =
N
E0 X
2
2
+
Ek cos
k(x a) + Fk sin
k(x a)
2
L
L
k=1
Teorema 1. Sea f una funcion continua a trozos en [a, b]. Entonces la suma parcial de la serie
de Fourier de f
N
2
2
a0 X
+
ak cos
k(x a) + bk sin
k(x a)
SN (x) =
2
L
L
k=1
da la aproximacion optima en el sentido de mnimos cuadrados de f (x) de entre todos los polinomios trigonometricos de grado a lo sumo N .
Esto significa que la diferencia
Z b
(f (x) TN (x))2 dx,
E=
a
con TN recorriendo los polinomios trigonometricos de grado menor o igual que N en [a, b] es
mnima cuando TN se elige la suma parcial SN (x) de la serie de Fourier de f .
La aproximacion se obtiene en el sentido de mnimos cuadrados; esto es, el polinomio PN era
la proyeccion ortogonal de f sobre el subespacio generado por las funciones
2
2
2
2
{1, cos
(x a) , sin
(x a) , . . . , cos
N (x a) , sin
N (x a) }, (3.33)
L
L
L
L
con el producto interno
Z
hf, gi =
f (x)g(x)dx.
a
Se puede comprobar que las funciones del conjunto (3.33) son ortogonales para este producto
interno, lo que lleva a las formulas (3.26) y (3.27) para los coeficientes. Esta idea de aproximaci
on
mnimos cuadrados motiva la aparicion de la serie infinita.
Teorema de convergencia
El teorema de convergencia fundamental para series de Fourier es el siguiente:
Teorema 2. Si f y f 0 son continuas a trozos en [a, b], entonces, para x (a, b), se tiene:
(1) Si f es continua en x,
a0 X
2
2
+
ak cos
k(x a) + bk sin
k(x a))
= f (x).
2
L
L
k=1
107
(2) Si f no es continua en x,
a0 X
2
2
f (x+) + f (x)
+
ak cos
k(x a) + bk sin
k(x a)
=
,
2
L
L
2
k=1
donde f(x+) y f(x) son, respectivamente, los lmites por la derecha y por la izquierda de
f en x.
Para x = a o x = b, la serie converge a 21 (f (a+) + f (b)).
Es decir, cuando f y f 0 son continuas trozos en [a, b], la serie de Fourier converge a f (x)
cuando f es continua en x y converge a la media entre los lmites a izquierda y derecha del punto
cuando f no es continua en x.
Hay que observar que la serie define una funcion periodica de perodo L.. Esto significa que
si extendemos f (x) del intervalo (a, b) a toda R por periodicidad, entonces la conclusion del
teorema se tiene para todo x sustituyendo f por su extension L-periodica.
Ejemplos.
(1) f : [, ] R,
X
1
f (x) = x 2
sin(k(x + )).
k
k=1
X
(1)k+1
k=1
sin(k(x + )).
1 2X 1
Seg
un el Teorema 2, si x (, 0)
1 2X 1
0=
sin((2k 1)(x + )).
2
2k 1
k=1
Y si x (0, ),
1 2X 1
1=
sin((2k 1)(x + )).
2
2k 1
k=1
3.6.3.
Diferenciaci
on e integraci
on de una serie de Fourier
Bajo ciertas condiciones, la serie de Fourier de una funcion puede derivarse termino a termino.
Teorema 3. Sea f continua en R y L-periodica. Supongamos que f 0 (x) y f 00 (x) existen y son
continuas a trozos en [a, b]. Entonces, la serie de Fourier de f (x) puede derivarse termino a termino
para producir la serie de Fourier de la derivada
2 X
2
2
h(x) = f (x)
k ak sin
k(x a) + bk cos
k(x a)
L
L
L
0
k=1
y si
x
Z
g(x) =
f (y)dy,
a
entonces
g(x) (x a)
a0 X
+
2
k=1
x
ak cos
2
2
k(y a) + bk sin
k(y a)
dy.
L
L
f (x) = x.
f (x) = x 2
g(x) =
x2
h(x) = 1 2
X
1
sin(k(x + )),
k
k=1
k=1
X
k=1
109
1
(1 cos(k(x + ))),
k2
cos(k(x + ))
Referencias
[1] J. W. Cooley, J W. Tukey, An algorithm for the machine calculation of complex Fourier
series, Math. Comput. 19, 297301 (1965).
[2] http://www.fftw.org.
[3] J. H. Mathews, K. D. Fink, Metodos numericos con MATLAB, Prentice Hall, 2000.
[4] P. J. Olver, Applied Mathematics Lecture Notes, http://www.math.umn.edu/olver/appl.html.
[5] J. M. Sanz Serna, Diez Lecciones de Calculo Numerico, Universidad de Valladolid, 1998.
110
Tema 4
M
etodos iterativos
1. Introduccion.
2. Metodos iterativos para sistemas lineales.
3. Metodos iterativos para ecuaciones no lineales.
4. Metodos iterativos para sistemas no lineales.
5. Ejercicios.
4.1.
Introducci
on
Esta leccion analiza la resolucion numerica de ecuaciones desde otro punto de vista. La primera
parte retoma el planteamiento del tema 2. Como ya se comento all, la b
usqueda de metodos de
resolucion de sistemas lineales grandes y con estructura (procedentes de aproximaciones por mnimos cuadrados) llevo a la construccion de metodos iterativos como alternativa a los directos. Tales
metodos son aproximados: generan una sucesion de vectores que, en el mejor de los casos, converge a la solucion del sistema. Las ventajas que presentan son, sin embargo, varias: en terminos
computacionales, se adaptan mucho mejor a la estructura del problema que los directos. Adem
as,
a pesar de ser aproximados, en muchas situaciones basta con un n
umero relativamente bajo de
aproximaciones para obtener una solucion bastante precisa. El aprovechamiento de las ventajas
de dispersion y estructura de las matrices que proporcionan estos metodos hace que su utilizaci
on
sea muy frecuente en la resolucion de sistemas lineales en circuitos electricos o en otros sistemas
de se
nales. Siguiendo un orden cronologico, presentaremos la formulacion de las tecnicas iterativas clasicas, analizando su implementacion. Por u
ltimo, al tratarse de metodos aproximados, es
obligado realizar un estudio sobre las condiciones a exigir para su convergencia a la soluci
on del
sistema.
La segunda parte de la leccion analiza el caso de tener ecuaciones no lineales, esto es, el
problema de aproximar ceros de funciones. La necesidad de obtener metodos de aproximaci
on
queda evidente con el simple comentario de que para un polinomio de grado cinco no hay formula
que calcule analticamente sus races. De manera que en general no pueden obtenerse exactamente
races de ecuaciones.
111
Desde el punto de vista de las aplicaciones, muchos problemas pueden modelarse matem
aticamente como una ecuacion, generalmente no lineal
f (x) = 0,
(4.1)
con f cierta funcion o sistema de funciones de una variable o conjunto de variables x, y cuya
resolucion analtica no es posible, o es posible pero demasiado elaborada.
Puede leerse una introduccion historica sobre la resolucion de ecuaciones por aproximaciones
sucesivas en [1]. Un ejemplo servira como ilustracion. El astronomo y matematico aleman Johannes Kepler (1571-1630), en uno de sus libros de Astronoma, necesitaba resolver la siguiente
ecuacion,
x = t e sin x,
con t y e valores dados. La ecuacion apareca al describir el movimiento de un cuerpo C que se
mueve en torno a otro A. De la primera ley de Kepler, C describe una elipse con A en uno de sus
focos. La variable e representa la excentricidad de la elipse, mientras que t juega el papel de un
angulo ficticio, que corresponde al movimiento de un cuerpo bajo movimiento circular uniforme.
Para obtener nuevas tablas astronomicas, Kepler necesitaba resolver la ecuacion para diferentes
valores de t. Llamo a su metodo la regla de posicion , en alusion probablemente a la larga
tradicion de las reglas de falsa posicion. El lo explica por medio de un ejemplo, con t = 50o 90 1000
y e = 1191000 . Razona como sigue: puesto que x < t, entonces Sin (x) < Sin (t) = 76775 (el
usa Sin para referirse al seno multiplicado por el radio 100000). Kepler elige un valor de x1
tal que Sin (x1 ) = 70000 con lo que obtiene x1 = 44o 250 (o mas precisamente 44o 250 3700 ). A
continuacion calcula
x1 + e sin x1 46o 440 ,
de donde
t (x1 + e sin x1 ) 3o 250 .
Ahora usa este error para obtener una mejor aproximacion a x:
x2 = x1 + 3o 250 = 47o 500 .
Puesto que Sin (x2 ) es cercano a 74000, hace los calculos mas sencillos tomando
x2 = ArcSin (74000) 47o 440 600 ,
(el valor correcto es 47o 430 5300 ). Usando este nuevo valor de x2 , Kepler repite los calculos,
x2 + e sin x2 50o 100 5900 ,
t (x1 + e sin x1 ) 1o 490 ,
obteniendo una nueva aproximacion
x3 = x2 1o 490 = 47o 420 1700 .
y parando aqu el algoritmo, puesto que x3 + e sin x3 = 50o 90 700 . El metodo es conocido hoy
como la iteracion de punto fijo.
En los Principia (1687), Newton mostro interes en la ecuacion de Kepler. En su libro, demostro como encontrar la posicion de un cuerpo moviendose en una elipse dada en un tiempo
112
asignado a traves de una prueba geometrica, pero reconocio que la construccion era difcil, proponiendo el uso de una solucion aproximada. Esta descripcion fue la primera publicaci
on de
Newton de lo que hoy llamamos metodo de Newton.
La idea de las tecnicas numericas para resolver este problema consiste en, dada una aproximacion inicial a un cero de la funcion, generar una sucesion de aproximaciones de forma algortmica. El proceso se detiene cuando se considera que el u
ltimo termino calculado de la sucesi
on es
una aproximacion suficientemente buena, en un sentido determinado por el algoritmo.
A manera de introduccion, discutiremos primero los metodos mas sencillos y menos eficientes,
para luego estudiar los metodos iterativos de punto fijo y de Newton, que son los mas utilizados
en la practica.
4.2.
M
etodos iterativos para sistemas lineales
4.2.1.
A Mnn (C), x Cn , b Cn ,
(4.2)
T
ecnicas iterativas cl
asicas
Los metodos iterativos generan, a partir de una aproximacion inicial a la solucion, una sucesi
on
de vectores que, en ausencia de errores de redondeo, converge a la solucion del sistema. En la
practica, se calcula un n
umero finito de iterantes (estableciendo un control para la iteraci
on) y
el u
ltimo calculado se toma como aproximacion a la solucion.
Las tecnicas iterativas que aqu mencionaremos parten de considerar una matriz Q, llamada
de descomposicion, escribir A = Q (Q A) y el sistema (4.2) en la forma
Qx = (Q A)x + b.
(4.3)
La ecuacion (4.3) genera una formula matricial recurrente. Partiendo de un dato inicial x(0)
arbitrario, se define la sucesion {x() } mediante la recurrencia
Qx() = (Q A)x(1) + b,
= 1, 2, . . . .
(4.4)
De manera que, conocido x(1) , el nuevo termino de la sucesion x() es la solucion del sistema
(4.4) con matriz Q y termino independiente (Q A)x(1) + b. Naturalmente, para que la tecnica
sea eficaz, es imprescindible que la eleccion de Q proporcione unos sistemas (4.4) mucho m
as
sencillos y baratos de resolver que el sistema original (4.2). Tambien, Q debe elegirse de modo
que la sucesion {x() } converja lo mas rapidamente posible a la solucion del sistema. Diferentes
elecciones de Q proporcionan las tecnicas iterativas clasicas que ahora comentaremos. Veremos
ademas que esta descomposicion permitira dar una recurrencia para estudiar los errores cometidos
y un criterio general de convergencia de los metodos.
4.2.2.
M
etodo de Jacobi
El metodo de Jacobi requiere que la matriz A del sistema (4.2) tenga elementos diagonales
no nulos, puesto que se toma como Q la matriz diagonal determinada por su diagonal principal:
Q = diag(a11 , . . . , ann ). De este modo, la formula (4.4) se escribe, por componentes,
X
()
(1)
aii xi = bi
aij xj
, i = 1, . . . , n,
j6=i
113
()
()
114
4.2.3.
M
etodo de Gauss-Seidel
Una modificacion del metodo de Jacobi viene dada al tomar como matriz Q de descomposici
on
la parte triangular inferior de A, incluyendo la diagonal
a11
a21
.
.
Q=
.
a
n1,1
an,1
0
a22
..
.
0
0
..
.
an1,2
an,2
an1,3
an,3
0
0
..
.
an1,n1
an,n1
0
0
..
.
,
0
ann
de manera que, a diferencia del metodo de Jacobi, se van incorporando a la iteracion las componentes del iterante nuevo inmediatamente despues de ser calculadas. Por componentes
()
aii xi
= bi
X
j<i
()
aij xj
(1)
aij xj
j>i
115
i = 1, . . . , n,
116
4.2.4.
M
etodo SOR
0
0
0
0
0
a21
0
0
0
0
.
.
.
.
..
..
..
..
E = ..
.
, F = D E A.
a
0
0
n1,1 an1,2 an1,3
an,1
an,2
an,3
an,n1 0
Entonces, el metodo de Gauss-Seidel se escribe, en forma matricial, como
Dx() = b + Ex() + F x(1) .
Partiendo de este metodo, definimos primero las componentes de un vector auxiliar x
()
()
aii x
i
= bi
i1
X
()
aij xj
j=1
n
X
(1)
aij xj
i = 1, . . . , n.
(4.5)
j=i+1
xi
(1)
= xi
()
+ (
xi
(1)
xi
(1)
) = (1 )xi
117
()
+
xi ,
i = 1, . . . , n.
(4.6)
i1
n
X
X
()
(1)
()
(1)
(1)
aii xi = aii xi
+ bi
aij xj
aij xj
aii xi
, i = 1, . . . , n.
j=1
j=i+1
O, matricialmente
(D E)x() = ((1 )D + F )x(1) + b,
es decir, que en la representacion (4.3), la matriz de descomposicion es
Q=
Al igual que la
un u
nico vector
La aparicion de
convergencia de
4.2.5.
1
(D E).
An
alisis de convergencia
Resultado general
Si escribimos (4.4) en la forma
x() = Q1 (Q A)x(1) + Q1 b,
y definimos M = Q1 (Q A), entonces los errores e = x() x satisfacen la recurrencia
e() = M e(1) ,
= 1, 2, . . . .
(4.7)
= D1 (E + F )
1
= (D E)
= (D E)
(Jacobi)
(Gauss-Seidel)
((1 )D + F ). (SOR)
(4.8)
Demostraci
on del teorema 4.2.1. Bastara ver que la sucesion de errores {e() }, que verifica
(4.7), converge a cero para cualquier vector inicial e(0) si y solo si (M ) < 1.
Supongamos primero que (M ) < 1. Sea un autovalor de M , con multiplicidad algebraica
m() y sea v un autovector generalizado asociado a . Recordemos la formula
m()1
M v=
X
k=0
k (M I)k v,
m()1
||M v||
X
k=0
(M )k ||wk ||.
||M v||
n
X
j vj ,
j=1
n
X
|j |||M vj || 0,
j=1
n
X
|aij |,
1 i n.
j=1,j6=i
Esto es, cada elemento de la diagonal principal es, en modulo, mayor que la suma de los modulos
de los restantes elementos de su fila. Por otro lado, una matriz real A es definida positiva si para
todo v 6= 0 v T Av > 0 o bien si todos sus autovalores son positivos.
119
m
.
log (M )
Cuando M es convergente ((M ) < 1), la cantidad R (M ) = log (M ) > 0 se llama velocidad
asintotica de convergencia y controla el n
umero de iteraciones necesario para reducir el error
inicial en una cantidad dada si dicho error es un autovector de M asociado a un autovalor de
modulo maximo. Cuando e(0) es general, se tiene la cota procedente de (4.7)
||e() || ||M ||||e(0) ||,
de donde
||e() ||
||e(0) ||
!1/
||M ||1/ .
(4.9)
||
se llama velocidad promedio de convergencia en
Si ||M || < 1, el n
umero R(M ) = ln ||M
120
Estudio comparativo
Varias de las ventajas de los metodos iterativos frente a los directos suelen observarse mejor
en la resolucion de sistemas lineales grandes, dispersos y con estructura. Un ejemplo clasico viene
dado por la matriz J J
4 1
0 0
1
1 4 1 0
1
4
1
0
.
..
..
..
.. ..
.
.
.
.
.
A= .
.
.
.
..
..
.. ..
..
..
.
.
.
.
.
0 0
1 4 1
1
12
11
10
N
s
7
s
S
Estimamos el potencial en un n
umero finito de puntos del aislante. Superponemos, como
ilustracion elemental, una malla a la seccion de la placa y buscamos el potencial en cada uno de
sus nudos, numerados en la figura desde el 1 hasta el 12. En la practica, la malla constar
a de
bastantes mas puntos.
La hipotesis fundamental del modelo es que en el estado estacionario, el potencial en cada
nudo de la malla es la media de los potenciales de los cuatro nudos mas proximos, esto es
1
VC = (VN + VS + VE + VO ).
4
La aplicacion de esta formula lleva a un sistema la matriz A anterior (para J = 12). El termino
independiente contiene valores del potencial en los puntos de la malla que caen en las placas.
Recordemos lo mencionado en secciones previas, en referencia a la implementacion. Un programa que implementa el metodo de Jacobi consta esencialmente de un lazo en el que se calcula
el iterante n-esimo a partir del anterior. A este lazo precedera la entrada del segundo miembro
y del iterante inicial. Como ya se ha comentado previamente, en cada momento, solo se guardan
en memoria dos iterantes: el que se esta calculando, que se guarda en un vector xn, y el anterior,
121
almacenado en, por ejemplo, xv. Al final del ciclo, el contenido de xn se sobreescribe en xv,
actualizandose para una nueva iteracion.
As, las componentes de xn verifican, salvo la primera y la u
ltima,
4 xn(i) = b(i) (xv(i 1) + xv( i + 1)).
Las ecuaciones primera y u
ltima son
4 xn(1) = b(1) (xv(2) + xv(J))
4 xn(J) = b(J) (xv(1) + xv(J 1)).
Sobreescribiendo 0,25 b en b, tenemos el ciclo
xn(1) = b(1) + 0,25 (xv(2) + xv(J))
i = 2:J 1
xn(i) = b(i) + 0,25 (xv(i 1) + xv(i + 1))
xn(J) = b(J) + 0,25 (xv(1) + xv(J 1)).
El lazo principal se ejecuta, por ejemplo, hasta que dos iterantes consecutivos difieran, en norma del maximo, en menos de una tolerancia, o bien dentro de otro ciclo que determina el n
umero
maximo de iteraciones. Hay que observar en este ejemplo la importancia de la estructura dispersa
de la matriz A. De este modo, el programa resuelve el sistema Ax = b, donde puede comprobarse que la matriz A es simetrica y definida positiva. Ademas, es diagonalmente dominante.
De manera que los metodos de Jacobi y Gauss-Seidel, para este sistema, son convergentes.
En cuanto al metodo de Gauss-Seidel, aqu, a diferencia del anterior, solo es necesario, dada
la construccion del algoritmo, un vector x, que esta guardando a la vez el iterante anterior y el
actual:
x(1) = b(1) + 0,25 (x(2) + x(J))
i = 2:J 1
x(i) = b(i) + 0,25 (x(i 1) + x(i + 1))
x(J) = b(J) + 0,25 (x(1) + x(J 1)).
Como ilustracion, tomemos un termino independiente b con todas sus componentes iguales a
2. En ese caso, la solucion del sistema es el vector con todas sus componentes iguales a uno.
Para J = 100 y distintas tolerancias, la ejecucion de los dos metodos da lugar a los siguientes
resultados:
M
etodo de Jacobi
122
T OL
N IT ER
||U u||
1E-05
1E-07
1E-09
1E-11
EPS
17
24
30
37
52
7.6294E-06
5.9605E-08
4.3132E-10
7.2759E-12
EPS
M
etodo de Gauss-Seidel
T OL
N IT ER
||U u||
1E-05
1E-07
1E-09
1E-11
EPS
12
16
20
24
34
1.8992E-06
2.3291E-08
2.8701E-10
3.5415E-12
EPS
en sistemas dispersos, siendo en general preferible el uso de metodos iterativos. Hay que indicar,
ademas, que la grafica muestra que el n
umero de operaciones del algoritmo tridiagonal y de los
metodos iterativos crece linealmente con el tama
no del sistema, mientras que el algoritmo general
aporta un n
umero de operaciones que crece con J 3 . Esto se observa estimando las pendientes de
las rectas asociadas a los metodos, teniendo en cuenta que la escala es logartmica. Los dos metodos iterativos tambien generan un n
umero de operaciones que crece linealmente con el tama
no
del sistema.
4.3.
4.3.1.
M
etodos iterativos para ecuaciones no lineales
M
etodo de bisecci
on
124
Algoritmo 4.5. M
etodo de bisecci
on.
Datos de entrada:
a, b: Extremos del intervalo.
Se supone que f (a)f (b) < 0. Puede incorporarse
un mensaje de error por si esto no ocurre.
tol: Tolerancia del proceso.
maxit: N
umero maximo de iteraciones.
niter = 0, error = 1
Mientras niter < maxit&error > tol
p1 (b a)/2
p a + p1
Evaluar fp f (p), fa f (a)
Si fp = 0
finalizar
error = |p1 |
niter niter + 1
Si fa fp > 0
ap
en otro caso
bp
Ejemplo 4.1. Tomando la ecuacion de Kepler mencionada en la introduccion, con los valores
t = 0,6, e = 0,5, la tabla 4.1 muestra los resultados del algoritmo de biseccion para los primeros
quince iterantes. La columna de errores se ha construido comparando el correspondiente iterante
con el u
ltimo calculado. La figura 4.2 es una ilustracion grafica del procedimiento.
iteracion
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
iterantes
0.7854
1.1781
0.9817
1.0799
1.0308
1.0063
1.0186
1.0247
1.0277
1.0293
1.0285
1.0282
1.0283
1.0282
1.0282
error
2.4280E-01
1.4990E-01
4.6452E-02
5.1722E-02
2.6351E-03
2.1909E-02
9.6368E-03
3.5010E-03
4.3287E-04
1.1011E-03
3.3412E-04
4.9377E-05
1.4237E-04
4.6497E-05
1.4400E-06
f(iterante)
-0.1682
0.1162
-0.0340
0.0390
0.0020
-0.0161
-0.0071
-0.0026
-0.0003
0.0008
0.0003
0.0000
0.0001
0.0000
0.0000
125
Hay dos comentarios a realizar sobre el algoritmo. El primero es que funciona siempre, aunque
la convergencia del procedimiento suele ser muy lenta, en comparacion con otros metodos. Otro
aspecto es que proporciona una estimacion directa del error cometido, que viene dada por el
tama
no del u
ltimo intervalo calculado.
Teorema 4.3.2 (Convergencia del metodo de bisecci
on). Sea f C([a, b]) con f (a)f (b) < 0.
ba
,
2n+1
n = 0, 1, . . .
Demostraci
on. Puede leerse en [2].
4.3.2.
M
etodo de la secante
126
conocidas aproximaciones xn , xn+1 , la recta que interpola a f en esos puntos tiene por ecuaci
on
y = f (xn+1 ) +
f (xn+1 ) f (xn )
(x xn+1 ).
xn+1 xn
f (xn+1 )
,
f [xn , xn+1 ]
f [xn , xn+1 ] =
f (xn+1 ) f (xn )
.
xn+1 xn
(4.10)
A pesar de que la convergencia no esta asegurada, cuando esta se da, el metodo es mas r
apido
que el de biseccion.
Ejemplo 4.2. Para el mismo problema del ejemplo 4.1, la tabla 4.2 muestra una mayor velocidad en la convergencia. La figura 4.3 es una ilustracion grafica del procedimiento.
iteracion
1
2
3
4
5
6
iterantes
0.9666
1.0327
1.0281
1.0282
1.0282
1.0282
error
0.0616
0.0045
0.0001
5.2483e-05
2.1301e-05
7.4054e-06
Algoritmo 4.6. M
etodo de la secante.
Datos de entrada:
p0 , p1 : Aproximaciones iniciales.
tol: Tolerancia del proceso.
maxit: N
umero maximo de iteraciones.
niter = 0, error = 1
Mientras niter < maxit y error > tol
Evaluar f0 f (p0 ), f1 f (p1 )
p p1 f1 (p1 p0 )/(f1 f0 )
error = |p p1 |
k k+1
p 0 p1
p1 p
El metodo de la secante dado por (4.10) puede completarse con un proceso de decision para
asegurar la convergencia [2]. Conocidos (xn , f (xn )) y (xn+1 , f (xn+1 )), recordemos que xn+2 es la
componente horizontal del punto de corte de la recta secante que pasa por tales puntos con el
eje horizontal. Supongamos f (xn )f (xn+1 ) < 0. Calculado (xn+2 , f (xn+2 )), se tienen tres posibilidades:
127
1. f (xn )f (xn+2 ) < 0, en cuyo caso hay un cero en [xn , xn+2 ]. Se eligen los puntos xn , xn+2
para continuar el proceso.
2. f (xn+1 )f (xn+2 ) < 0, en cuyo caso hay un cero en [xn+1 , xn+2 ]. Ahora se eligen los puntos
xn+1 , xn+2 para continuar el proceso.
3. f (xn+2 ) = 0. No hay nada mas que hacer, porque ya hemos encontrado la raz.
Esta variante genera un nuevo metodo, conocido como regula falsi o falsa posicion. Adem
as,
permite construir una sucesion de intervalos [an , bn ] que contienen un cero de f . La aproximaci
on
viene dada por
cn = bn
f (bn )
,
f [bn , an ]
4.3.3.
M
etodo de punto fijo
Descripci
on
Dos metodos de aproximacion mas eficientes son la iteracion de punto fijo y el metodo de
Newton.
Ya hemos comentado algo del primero en la introduccion. La iteracion de punto fijo se basa
en reescribir el problema (4.1) de hallar un cero de f en el problema de calcular un punto fijo de
128
(4.11)
n = 1, 2, . . . .
(4.12)
podemos asegurar la convergencia de (4.12). Esta
exige, ademas, condiciones sobre la funci
on g
que influyen en su eleccion. Todo ello se recoge en el siguiente resultado.
Teorema 4.3.4 Supongamos que g es de clase C 1 en un intervalo (r , r + ), con > 0. Si
|g 0 (r)| < 1,
entonces, existe 0 < < tal que si x0 [r , r + ] se puede construir la sucesi
on (4.12) con
xn [r , r + ] y convergente a r. Adem
as, si los errores en = xn r son no nulos, los cocientes
verifican
en+1
g 0 (r),
en
n .
Nota: Dado que no conocemos r, el teorema es un resultado teorico que establece las condiciones
para la convergencia de (4.12), aunque en la practica no puede utilizarse.
Demostraci
on. Sea M con |g 0 (r)| < M < 1. Como g 0 es continua, existe < tal que
|g 0 (x)| M,
x [r , r + ].
(4.13)
Por tanto
|xn+1 r| = |en+1 | = |g 0 (yn )en | M |xn r| < |xn r| .
Luego xn+1 [r , r + ]. Ademas, de (4.13),
|en+1 | M |en | M 2 |en1 | M n |e0 |,
con M < 1, lo que prueba que en 0 y por tanto la sucesion {xn }n converge a r. Por u
ltimo,
tambien de (4.13), si en 6= 0 y n , como xn r, tambien yn r. Al ser g 0 continua
en+1
= g 0 (yn ) g 0 (r).
en
Ejemplo 4.3. Para el mismo problema del ejemplo 4.1, la tabla 4.3 muestra, en su segunda
columna, la sucesion de iterantes obtenidos mediante punto fijo. Con cuatro cifras decimales, la
iteracion converge hacia 1,0282. El cociente entre dos errores consecutivos (
ultima columna) se
aproxima tambien a una cantidad constante que, siguiendo el teorema 4.3.4 es g 0 (r) 0,2582.
130
n
2
3
4
5
6
7
8
9
10
11
12
13
xn
0.600
0.8823
0.9861
1.0169
1.0253
1.0274
1.0280
1.0281
1.0282
1.0282
1.0282
1.0282
en
4.2818E-01
1.4586E-01
4.2073E-02
1.1238E-02
2.9286E-03
7.5797E-04
1.9582E-04
5.0567E-05
1.3056E-05
3.3711E-06
8.7037E-07
2.2472E-07
en+1 /en
0.3407
0.2884
0.2671
0.2606
0.2588
0.2583
0.2582
0.2582
0.2582
0.2582
0.2582
Tabla 4.3: Iteracion de punto fijo para f (x) = x 0,5 sin x 0,6 = 0, con g(x) = 0,5 sin x + 0,6.
131
Algoritmo 4.7. M
etodo de punto fijo.
Datos de entrada:
p0 : Aproximacion inicial.
tol: Tolerancia del proceso.
maxit: N
umero maximo de iteraciones.
niter = 0, error = 1
Mientras niter < maxit y error > tol
Evaluar p1 g(p0 )
error = |p1 p0 |
k k+1
p 0 p1
El teorema 4.3.4 afirma, entre otras cosas, que cuando la iteracion de punto fijo es convergente,
su velocidad de convergencia es lineal. Vamos a precisar esto un poco. Sea {xn }n una sucesi
on
que converge a un punto x y sea en = xn x. Si existen A, B > 0 tales que
lm
|en+1 |
= A,
|en |B
entonces se dice que la sucesion converge a x con orden de convergencia dado por B. El n
umero A
se llama constante asintotica del error. En particular, el caso B = 1 se llama convergencia lineal
y el caso B = 2, convergencia cuadratica.
De manera que el teorema 4.3.4 dice que, en caso de convergencia de la iteracion de punto
fijo, esta es tpicamente lineal con constante asintotica g 0 (r). Un caso a
un mas favorable puede
darse en la siguiente situacion.
Teorema 4.3.5 Supongamos que, en las condiciones del teorema 4.3.4, g 0 (r) = 0 y que existe
g 00 (r). Entonces, la iteraci
on de punto fijo converge con orden cuadr
atico y constante asint
otica
g 00 (r)/2.
Demostraci
on. Si en la demostracion del teorema 4.3.4 desarrollamos por Taylor
en+1 = xn+1 r = g(xn ) g(r) = g 0 (r)(xn r) +
g 00 (r)
(xn r)2 + o((xn r)2 ),
2
entonces, si en = xn r 6= 0,
en+1
g 00 (r) o(e2n )
g 00 (r)
=
+
.
e2n
2
e2n
2
El estudio local de la convergencia se completa con los casos que dan lugar a divergencia local.
Teorema 4.3.6 Supongamos que g es de clase C 1 en un intervalo (r , r + ), > 0, en torno
a un punto fijo de g. Si
|g 0 (r)| > 1,
entonces, las u
nicas sucesiones de la forma (4.12) convergentes a r son aquellas tales que xn = r
a partir de un n.
132
Demostraci
on. Sea M con |g 0 (r)| > M > 1. Como g 0 es continua, existe < tal que
|g 0 (x)| M,
x [r , r + ].
4.3.4.
El m
etodo de Newton
Descripci
on
Ideado por Newton a mediados del siglo XVII, es tambien conocido como metodo de la
tangente [1]. Fue formulado en los terminos en los que se conoce actualmente por Raphson en 1690,
incluyendo su construccion geometrica. Supongamos que la aproximacion inicial x0 esta cerca de
la raz r y supongamos que f es derivable. Como segunda aproximacion x1 se toma la intersecci
on
del eje de abscisas con la recta tangente a la curva en (x0 , f (x0 )). La expresion de tal recta es
y f (x0 ) = f 0 (x0 )(x x0 ) de modo que
x1 = x0
f (x0 )
,
f 0 (x0 )
f (xn )
,
f 0 (xn )
n = 0, 1, . . . .
(4.14)
La relacion con la iteracion de punto fijo se encuentra en la interpretacion de (4.14) como dicha
iteracion para la funcion
g(x) = x
f (x)
.
f 0 (x)
Esta conexion nos permite entonces deducir aspectos de convergencia del nuevo metodo.
133
Algoritmo 4.8. M
etodo de Newton.
Datos de entrada:
p0 : Aproximacion inicial.
tol: Tolerancia del proceso.
maxit: N
umero maximo de iteraciones.
niter = 0
Valor inicial de error
Mientras niter < maxit y error > tol
Evaluar f0 = f (p0 ), df0 = f 0 (p0 )
p1 p0 f0 /df0
error = |p1 p0 |
k k+1
p 0 p1
4.3.5.
Aspectos de convergencia
f (x)f 00 (x)
,
f 0 (x)2
134
de manera que, con las hipotesis del teorema 4.3.7, se tiene que g 0 (r) = 0. Se deduce entonces el
resultado usando el teorema 4.3.5, que indica ademas que la constante asintotica es g 00 (r)/2.
Al igual que en el caso de la iteracion de punto fijo, el metodo de Newton requiere tambien
una buena aproximacion inicial al cero r para la convergencia. Sin embargo, cuando esta se da,
es tpicamente mas rapida. El metodo tambien exige mayor regularidad de la funcion.
La condicion f (r) = 0, f 0 (r) 6= 0 del teorema 4.3.7 significa que r ha de ser un cero simple de
f . Para ceros m
ultiples, la convergencia del metodo ya no es cuadratica (veanse los ejercicios al
final del tema).
Ejemplo 4.4. Para el mismo problema del ejemplo 4.1, la tabla 4.4 muestra, en su segunda
columna, la sucesion de iterantes obtenidos mediante el metodo de Newton. Ahora los cocientes
en+1 /e2n (
ultima columna) muestran la convergencia cuadratica. El n
umero de iteraciones es
tambien bastante reducido, en comparacion con los otros procedimientos.
n
0
1
2
3
xn
1.2000
1.0364
1.0282
1.0282
en
1.7182E-01
8.1936E-03
1.9347E-05
1.0803E-10
en+1 /e2n
0.2775
0.2882
0.2886
135
4.4.
M
etodos iterativos para sistemas no lineales
(4.15)
fm (x1 , . . . , xm ) = 0.
Para describir tanto la iteracion de punto fijo para sistemas como, en la seccion siguiente, el
metodo de Newton para sistemas, necesitamos manejar normas matriciales y sus propiedades,
explicadas en el tema 2, dado que ahora el papel de la derivada vendra representado por la
matriz jacobiana.
4.4.1.
Iteraci
on de punto fijo para sistemas
n = 1, 2, . . .
Denotamos por || || una norma en Rm . Los siguientes resultados establecen las condiciones
para la convergencia en el caso de sistemas. Denotamos por g = (g1 , . . . , gm ) a la funci
on de
iteracion que tenga por punto fijo al cero r Rm de f a aproximar. En estas condiciones, g 0 (r)
es la matriz jacobiana de g en el punto r.
Teorema 4.4.1 Supongamos que g est
a definida y es de clase C 1 en un entorno de r, con
0
||g (r)|| < 1, donde || || tambien representa la norma matricial asociada a la vectorial definida en
Rm . Entonces, existe tal que si ||x0 r|| < , la sucesi
on xn = g(xn1 ) verifica que ||xn r|| <
para todo n y es convergente a r . Adem
as, si {xn } converge a r, existe 0 < K < 1 tal que los
errores en = xn r verifican
||en+1 || K||en ||,
||en+1 g 0 (r)en || = o(1)||en ||,
n .
4.4.2.
El m
etodo de Newton para sistemas
La deduccion del metodo de Newton para obtener un cero r = (r1 , . . . , rm ) del sistema (4.15),
donde suponemos que f = (f1 , . . . , fm ) es de clase C 1 en un disco abierto de Rm , es similar a
136
la del metodo para ecuaciones. Conocida la aproximacion xn , la siguiente xn+1 se obtiene como
interseccion de la variedad y = 0 y el hiperplano tangente a y = f (x) en (xn , f (xn )), de manera
que
0 = f (xn ) + f 0 (xn )(xn+1 xn ),
(4.16)
y, por tanto, xn+1 satisface un sistema lineal con matriz dada por f 0 (xn ). Para la existencia
de xn+1 se requiere entonces que la matriz jacobiana f 0 (xn ) sea invertible. La relacion con la
iteracion de punto fijo viene dada por la funcion
g(x) = x f 0 (x)1 f (x),
que tiene sentido cuando f 0 es invertible. La convergencia cuadratica se mantiene cuando f es de
clase C 3 en un entorno de r, f 0 (r) es invertible y la aproximacion inicial x0 se toma suficientemente
cerca de r [5].
Nota. Criterios de parada. En la programacion de los metodos que hemos mencionado en
esta leccion, pueden darse uno o varios criterios de parada del proceso. Los mas utilizados suelen
ser los siguientes (ya hemos hablado de alguno de ellos):
1. Establecer un n
umero maximo de iteraciones.
2. Fijar una tolerancia previa , suficientemente peque
na y parar el proceso en el paso n si
|f (xn )| < .
3. Detener las iteraciones cuando dos valores consecutivos xn y xn+1 estan suficientemente
cerca. Puede estimarse en terminos absolutos fijando una tolerancia T OL y parando el
proceso cuando
|xn+1 xn | < T OL,
o bien en terminos relativos, cuando
2|xn+1 xn |
< T OL.
|xn+1 | + |xn |
No es conveniente elegir tolerancias demasiado peque
nas. Una eleccion habitual consiste en
M
tomarlas unas cien veces mayores que 10 , con M el n
umero de cifras decimales de la
representacion en coma flotante que se utilice [2].
4.5.
Ejercicios
1 2 2
1 2
A1 = 1 1 1
,
A2 = 1 1
2 2 1
1 1
y b1 = b2 = (1, 2, 1)T .
137
1
0
1
xn+1 =
0
a
b
xn + f ,
(4.17)
donde es un parametro. Determina si es posible encontrar valores de con > a > b de forma
que el metodo (4.17) sea convergente.
Ejercicio 3. Estudia la convergencia del metodo
2 1 0
1 2 1
0 1 2 1
.
..
..
..
.
.
.
A = ..
.
.
.
.
.
.
.
.
.
0
0 0
0
0
0
0
0
0
..
..
.
.
.
..
..
..
.
.
.
1 2 1
0
aii = 1 >
|aij |,
1 i n.
j=1,j6=i
138
0 2 3
0
()
(1)
x = Mx
+ c, M =
1
0 0 , c = 0.
1 0 0
2
Estudia si la solucion dada por el metodo iterativo es la misma que la del sistema y si el metodo
es convergente. Determina cual de los cuatro metodos es, en principio, mas rapido.
Ejercicio 7. Demuestra que, si en denota el error en la iteracion n-esima del metodo de la secante,
entonces
en+1
f 00 (r)
= 0 .
n en en1
2f (r)
lm
Suponiendo que
|en+1 |
C,
|en |
deduce el orden del metodo.
Ejercicio 8. Comportamiento global de la iteraci
on de punto fijo [3]. Determina el
n
umero de puntos fijos de g(x) = (1/2)x x3 . Halla un punto > 0 con la propiedad g() = .
139
Estudia la iteracion de punto fijo para: (i) x0 (0, ). (ii) x0 = . (iii) x0 > . Que ocurre
cuando x0 es negativo?
Ejercicio 9. En el dibujo siguiente se muestra la grafica de cierta funcion g en el intervalo
[0, 5/2]. Los puntos A, B, C son los cortes de g con la bisectriz en dicho intervalo. Se sabe que
las pendientes de la grafica de g en los puntos A, B, C son, respectivamente, 1/4, 3/2 y 1/2.
Con estos datos, discute el comportamiento de la iteracion de punto fijo con la funcion g seg
un
la eleccion del iterante inicial.
Ejercicio 10. (i) Comprueba que, para [0, 4], la funcion g(x) = x(1 x) aplica el intervalo
[0, 1] en s mismo.
(ii) Sea x [0, 1], [0, 4]. Determina el n
umero de puntos fijos de g.
(iii) Comprueba que el punto fijo r = 0 es atractor si < 1 y repulsor si > 1. Comprueba que
para = 1 cualquier x0 [0, 1] conduce a una sucesion de iterantes que converge a r = 0 y que
cada error en es asintoticamente igual al precedente.
(iv) Estudia para que valores de el punto fijo no nulo es atractivo.
Ejercicio 11. Puntos fijos con |g 0 ()| = 1. Las funciones g(x) = sin x y g(x) = tan x tienen
ambas el punto fijo r = 0 y para ambas g 0 (0) = 1. Prueba que, para |x0 | suficientemente peque
no,
con el seno la iteracion de punto fijo converge mientras que con la tangente diverge. Por tanto el
caso |g 0 ()| = 1 es dudoso, la convergencia o divergencia depende de los valores de las derivadas
superiores de g.
Pruebe que si con |g 0 ()| = 1 hay convergencia al punto fijo entonces cada error es asint
oticamente de la misma magnitud del anterior, por lo que la convergencia es muy lenta y carece de
utilidad para hallar el punto fijo.
140
x2 + 2a
,
(1 + a)x
con a R un parametro, a 6= 1.. determina para que valores de a la iteracion de punto fijo
converge y para que valores la convergencia es cuadratica.
Ejercicio 13. (i) Si se utiliza el metodo de Newton con f (x) = x3 2, comenzando con x0 = 1.
Quien es x2 ?
(ii) Como podemos, utilizando el metodo de Newton-Raphson, implementar la divisi
on de
n
umeros reales? Y las races N -esimas de un n
umero?
Ejercicio 14. Considera una variacion del metodo de Newton en la que solamente se requiere
una derivada,
f (xn )
xn+1 = xn 0
.
f (x0 )
Analiza la convergencia de la sucesion de iterantes {xn } al cero x de f si se toma la sucesi
on de
iterantes suficientemente proxima al cero de f . Cual es la variacion equivalente del metodo de
Newton para sistemas?
Ejercicio 15. M
etodo de Newton para races m
ultiples. Sea f una funcion definida
un intervalo abierto que contenga un punto r con f (r) = 0. Supongamos que f es derivable
cualquier orden en todo punto y que r es el u
nico cero de f . Comprueba que si r es un cero
orden p > 1, entonces el metodo de Newton para f converge linealmente. Aplica el metodo
Newton para resolver x2 = 0 a partir de x0 dado y comprueba la convergencia lineal.
en
de
de
de
pf (xn )
,
f 0 (xn )
2 +y 2 +8
10
, xy
2 +x+8
10
Ejercicio 18. Determine una region en el plano tal que la iteracion de punto fijo aplicada al
sistema
x = g1 (x, y) = (x2 y 2 x 3)/3,
y = g2 (x, y) = (x + y + 1)/3,
sea convergente para cualquier punto inicial de dicha region.
Ejercicio 19. Para la resolucion de un sistema lineal Que relacion hay entre el metodo de Jacobi
y la iteracion de punto fijo aplicada al sistema?
141
Referencias
[1] J.-L.Chabert (Ed.), A History of Algorithms. From the Pebble to the Microchip, Springer,
1999.
[2] J. H. Mathews, K. D. Fink, Metodos numericos con MATLAB, Prentice Hall, 2000.
[3] J. M. Sanz Serna, Diez Lecciones de Calculo Numerico, Universidad de Valladolid, 1998.
[4] G. W. Stewart, Afternotes on Numerical Analysis, SIAM, 1998.
[5] G. W. Stewart, Afternotes Goes to Graduate School: Lectures on Advanced Numerical
Analysis, SIAM 1999.
142
Tema 5
Cuadratura y derivaci
on num
erica
1. Introduccion.
2. Metodos elementales de construccion de reglas de cuadratura. Analisis del error.
3. Reglas compuestas.
4. Formulas de aproximacion a la derivada.
5. Ejercicios.
5.1.
Introducci
on
El tema 3 abordaba diferentes maneras de aproximar funciones o de representar datos utilizando interpolacion. En esta leccion trataremos otro problema basico en analisis numerico, que
es la aproximacion a integrales y derivadas. La cuadratura numerica y la derivacion numerica
comparten argumentos para su justificacion, as como metodos de construccion. La cuadratura
numerica es el proceso por el que se genera un valor numerico para la integracion de una funcion sobre un conjunto. Muchas reglas proceden del calculo aproximado de areas y vol
umenes
realizado en tiempos antiguos [2]. La invencion del calculo infinitesimal a finales del siglo XVII
establecio la interpretacion de estos en terminos de integrales, de manera que las antiguas reglas
de aproximacion fueron generalizadas a traves de una b
usqueda de reglas de cuadratura para
integrales definidas. Hay ejemplos de estas que no pueden calcularse con las tecnicas elementales
de calculo (dependientes de la obtencion de primitivas): tpicamente, para encontrar el valor de
Z b
I(f ) =
f (x)dx,
(5.1)
a
uno debe obtener primero una funcion F con F 0 = f y aplicar la regla de Barrow: I(f ) =
F (b) F (a). Hay, sin embargo, funciones mas o menos elementales que carecen de primitivas
2
sencillas. Quiza el ejemplo mas utilizado es el de f (x) = ex . Una primitiva de esta funcion es [3]
F (x) =
X
k=0
x2k+1
,
(2k + 1)k!
143
que no es abordable numericamente. Una idea para la aproximacion numerica a (5.1) consiste en
sustituir f por otra funcion g que aproxime a f de manera adecuada y sea facil de integrar
b
f (x)dx
g(x)dx.
a
En este sentido, los polinomios son funciones que pueden satisfacer los dos requisitos, si bien
el tipo de aproximacion puede ser diferente, en funcion del problema de interpolacion de f que
resuelva el polinomio. Por ejemplo, g puede proceder de la truncacion de la serie de Taylor. Para
el caso anterior
Z
ex dx
N
1X
0 k=0
X
x2k
1
dx =
.
k!
(2k + 1)k!
k=0
1730
5.2.
M
etodos elementales de construcci
on de reglas de cuadratura
N
X
j f (xj ).
(5.2)
j=0
Supondremos a y b finitos, con f continua a trozos en [a, b]. En la combinacion (5.2) de evaluaciones del integrando, los puntos de evaluacion xj [a, b], j = 0, . . . , N se llaman nodos o abscisas
de cuadratura. Los coeficientes j de la combinacion son los llamados pesos o coeficientes de la
formula.
144
Punto kilometrico
1730
1735
1740
1745
1750
1755
1760
1765
1770
1775
1780
1785
1790
1795
1800
1805
1810
5.2.1.
M
etodo interpolatorio
Sea PN el u
nico polinomio de grado menor o igual que N que interpola a f en los nodos de
cuadratura x0 , . . . , xN . En forma de Lagrange, se escribe (tema 3)
PN (x) =
N
X
j=0
145
Z
f (x)dx
PN (x)dx =
a
N Z
X
j=0
Lj (x)dx f (xj ),
(5.3)
Lj (x)dx,
j = 0, . . . , N.
No es difcil ver que una formula interpolatoria (5.3) tiene un grado de precision mayor o igual
que N , pues si f es un polinomio de grado menor o igual que N , coincide con su polinomio
interpolador.
Interpolatorias son las cuatro reglas clasicas:
1. Regla del rectangulo: N = 0, x0 = a.
IREC (f ) = (b a)f (a).
Grado de precision cero.
2. Regla del punto medio: N = 0, x0 =
a+b
2 .
IP M (f ) = (b a)f
a+b
2
.
(b a)
(f (a) + f (b)) .
2
146
I(f )
I(f )
I(f )
I(f )
5.2.2.
(b a)
(f (a) + f (b))
2
(b a)
a+b
f (a) + 4f
+ f (b)
6
2
ba
a + 2(b a)
3(b a)
f (a) + 3f a +
+ 3f
+ f (b)
8
3
3
2(b a)
ba
a + 2(b a)
7f (a) + 32f a +
+ 12f
45
4
4
3(b a)
+32f a +
+ 7f (b)
4
M
etodo directo
El metodo consiste en, conocidos los nodos x0 , . . . , xN , determinar los coeficientes imponiendo
directamente las condiciones para obtener un grado de precision, haciendo la regla exacta para
una base de polinomios del grado elegido. Notemos que la regla (5.2) tiene grado de exactitud al
menos k si y solo si es exacta para f (x) = 1, x, . . . , xk . Esto lleva al sistema de k + 1 ecuaciones
para las N + 1 incognitas
0 + 1 + + N
= b a,
b2 a2
0 x0 + 1 x1 + + N xn =
,
2
..
..
.
.
bk+1 ak+1
0 xk0 + 1 xk1 + + N xkn =
.
k+1
(5.4)
b3 a3
,
3
que los anteriores valores de 0 , 1 no verifican, por lo que la regla obtenida tiene grado de
exactitud uno. Se trata de la regla de los trapecios antes presentada.
5.2.3.
M
etodo del desarrollo de Taylor
148
5.3.
An
alisis del error
(5.5)
que nos den una idea de la precision de la formula utilizada, as como un criterio de selecci
on de
las reglas.
La representacion de (5.5) que utilizaremos viene dada por el llamado Teorema de Peano y la
aplicaremos tambien en la seccion siguiente al tratar las reglas compuestas. La demostraci
on del
teorema del n
ucleo de Peano en su forma general puede consultarse en varias referencias [3, 1, 8].
Es asmismo aplicable en el analisis de errores de otros problemas, como por ejemplo la derivaci
on
numerica, que luego estudiaremos.
Teorema 5.3.1 Supongamos que la f
ormula (5.2) tiene grado de precisi
on al menos M . Entonces, para cada f C M +1 ([a, b]),
Z b
E(f ) =
f M +1) (t)K(t)dt,
a
donde K es la funci
on llamada n
ucleo de Peano, cuya expresi
on es
K(t) =
1
E(gt ),
M!
siendo gt la funci
on
gt (x) =
(x t)M
0
xt
x<t
Demostraci
on. Utilizando la formula de Taylor con resto integral
Z x
(x a)M M )
1
f (x) = f (a) + +
f (a) +
f M +1) (t)(x t)M dt,
M!
M! a
notese que, seg
un la definicion de gt , se puede escribir
Z x
Z b
f M +1) (t)(x t)M dt =
f M +1) (t)gt (x)dt.
a
Teniendo ahora en cuenta que la formula (5.2) tiene grado de precision al menos M
Z b
Z b
1
1
M +1)
M +1)
E
f
(t)gt (x)dt =
f
(t)E(gt )(x)dt .
E(f ) =
M!
M!
a
a
La u
ltima igualdad se debe a la forma de E, que le hace conmutar con la integral.
Corolario 5.3.2 Si en la situaci
on del teorema 5.3.1, el n
ucleo K no cambia de signo en [a, b],
entonces para cada f C M +1 ([a, b]) existe [a, b] tal que
E(f ) =
f M +1) ()
E(xM +1 ).
(M + 1)!
149
Demostraci
on. Del teorema del valor medio para el calculo integral [5] y el teorema 5.3.1,
se tiene
Z b
M +1)
K(t)dt,
E(f ) = f
()
a
E(x
Z
) = (M + 1)!
K(t)dt,
a
de donde
b
E(xM +1 )
.
(M + 1)!
K(t)dt =
a
Ejemplo 3.4. Error en la regla del rectangulo.
Z b
EREC (f ) =
f (x)dx (b a)f (a).
a
En este caso M = 0 y si f C 1 ,
b
Z
EREC (f ) =
f 0 (t)K(t)dt,
1
0
si
si
xt
x<t
Entonces
Z
Z
gt (x)dx (b a)gt (a) =
dx = b t,
t
y por tanto
Z
EREC (f ) =
f 0 (t)(b t)dt.
Ademas, como K(t) no cambia de signo en [a, b], existe [a, b] tal que
Z b
(b a)2 0
0
0
xdx (b a)a =
f ().
EREC (f ) = f ()EREC (x) = f ()
2
a
Ejemplo 3.5. Error en la regla del punto medio.
Z b
a+b
EP M (f ) =
f (x)dx (b a)f
.
2
a
150
f 00 (t)K(t)dt,
(x t)
0
si x t
si x < t
Entonces
Z
gt (x)dx (b a)gt
Z b
a+b
=
(x t)dx (b a)gt
2
t
(
2
K(t) = EP M (gt ) =
(bt)2
2
(b
(bt)
2
a) a+b
2
t =
(ta)2
2
a+b
2
si
tc=
si
a+b
2
tc
Observemos tambien que K(t) no cambia de signo en [a, b], luego existe [a, b] tal que
!
f 00 ()
f 00 () b3 a3
a+b 2
2
EREC (f ) =
EP M (x ) =
(b a)
2
2
3
2
=
5.4.
(b a)3 00
f ().
24
Reglas compuestas
La aparicion de las reglas de cuadratura compuestas se debe a que, como hemos visto en los
ejemplos anteriores, el error de cuadratura en las reglas simples depende de la longitud del intervalo de integracion. Para hacer el error peque
no, sera entonces necesario un intervalo peque
no.
Una alternativa para reducir el error consiste en dividir el intervalo [a, b] en una particion
z0 < z1 < < zN ,
aplicar la propiedad de aditividad de la integral
Z
f (x)dx =
a
N Z
X
i=1
zi
f (x)dx,
zi1
y aproximar cada una de las integrales por una regla de cuadratura (no tiene por que ser la misma
regla en cada subintervalo). As se obtienen las llamadas reglas de cuadratura compuestas. En
muchas ocasiones, la particion del intervalo se suele hacer uniforme, es decir que para todo
i = 0, . . . , N zi zi1 = h es constante.
Ejemplo 3.6. La aplicacion de una regla clasica en cada subintervalo de la particion proporciona la correspondiente regla compuesta. Por ejemplo:
151
5.4.1.
An
alisis del error en las reglas compuestas
En este caso, se puede analizar el error en las reglas compuestas utilizando el teorema 5.3.1
pero en cada subintervalo de la particion. Podemos ilustrar el procedimiento con la regla del
punto medio compuesta
Z b
N
X
zi1 + zi
f (x)dx
.
(zi zi1 )f
2
a
i=1
i=1
zi1
con f C 2 y
(
Ki (t) =
(zi t)2
2
(tzi1 )2
2
si
si
zi1 +zi
2
zi1 +zi
2
De este modo, si K(t) es la funcion real definida en [a, b] tal que K(t) = Ki (t) si t [zi1 , zi ],
entonces, si f C 2 [a, b]
Z b
E(f ) =
K(t)f 00 (t)dt.
a
As,
Z
|E(f )|
=
|f 00 (t)|dt
a
Z b
1
max |zi zi1 |2
|f 00 (t)|dt.
8 iiN
a
max |K(t)|
atb
(5.6)
Se dice entonces que el error tiene una convergencia de orden dos: la cota (5.6) sugiere que para
integrandos dos veces derivables con derivada segunda continua, la division por dos del diametro
h de la particion provoca que el error se divida aproximadamente por cuatro.
152
5.4.2.
Estudio comparativo
Para la eleccion de una regla de cuadratura en un problema concreto, hay varios criterios a
tener en cuenta. Es claro que uno de ellos ha de ser el de minimizar el error, de modo que es
adecuado utilizar reglas cuyo error de cuadratura tenga un alto grado de convergencia. Hay que
observar, en este sentido, que las cotas teoricas, como (5.6), no dependen solo del diametro h
de la particion; tambien es necesaria cierta regularidad del integrando. Sin este otro requisito, la
cota no es valida y no conocemos el comportamiento del error de cuadratura.
Otro criterio, a combinar con el anterior, tiene en cuenta el coste operativo. Reglas con
alto grado de convergencia pueden requerir muchas evaluaciones del integrando y, por tanto, ser
demasiado costosas computacionalmente. Hay entonces que establecer un cierto equilibrio entre
convergencia y coste operativo para la seleccion de la regla en un problema determinado.
En las siguientes tablas y figuras se ilustra una comparacion entre cuatro reglas compuestas,
rectangulo (RC), punto medio (PMC), trapecios (TC) y Simpson (SC), aplicadas en un intervalo
[0, 1], con una particion uniforme de longitud h. Vamos a ilustrar su comportamiento con dos
integrandos de diferente regularidad.
h
5.0e-001
2.5e-001
1.25e-001
6.25e-002
RC
2.1998e-001
1.0758e-001
5.3191e-002
2.6446e-002
PMC
4.8237e-003
1.1993e-003
2.9942e-004
7.4829e-005
TC
9.6172e-003
2.3968e-003
5.9872e-004
1.4965e-004
SC
1.0051e-005
6.2467e-007
3.8987e-008
2.4358e-009
con las cuatro reglas antes mencionadas. La figura 5.1 muestra el error (en escala doblemente
logartmica) en funcion del diametro h para cada regla. La figura 5.2 ilustra, tambien en escala
doblemente logartmica, el error de cuadratura frente a las evaluaciones del integrando realizadas
por cada regla, para los valores de h de la tabla 5.2.
h
5.0000e-001
2.5000e-001
1.2500e-001
6.2500e-002
RC
3.1311e-001
1.4838e-001
7.1036e-002
3.4335e-002
PMC
1.6346e-002
6.3107e-003
2.3655e-003
8.7001e-004
TC
SC
6.3113e-002 1.0140e-002
2.3384e-002 3.5874e-003
8.5364e-003 1.2685e-003
3.0855e-003 4.4848e-004
Podemos deducir alguna informacion de estos resultados. Notemos primero que el integrando
f (x) = sin x es una funcion derivable de cualquier orden. De los datos de la tabla 5.2 podemos hacernos una idea sobre la convergencia de los metodos. Por ejemplo, cuando dividimos el
diametro h entre dos (dos filas consecutivas de la tabla), los errores de (RC) se dividen, aproximadamente, por dos, los de (PMC) y (TC) por cuatro y los de (SC) por dieciseis. Esto sugiere
153
anteriores pero con el integrando f (x) = x, x [0, 1], que es solo continuo. Ahora, la tabla 5.3 y
la figura 5.3 muestran una perdida de convergencia de los metodos (PMC), (TC) y (SC), siendo
ahora comparables a (RC), si bien para un valor fijo de h, los errores menores siguen siendo
proporcionados por (SC). Por su parte, la figura 5.4 indica que, para los valores de la tabla 5.3, la
regla (PMC) se presenta ahora mas competitiva, aunque, en comparacion con el ejemplo anterior,
los errores obtenidos son mayores.
154
Figura 5.2: Comparacion de errores con evaluaciones del integrando para f (x) = sin x.
5.5.
F
ormulas de aproximaci
on a la derivada
Buscamos ahora formulas numericas para aproximar la derivada de una funcion f en un punto
y en la forma
DN +1 =
N
X
j f (xj ),
(5.7)
j=0
f (y)
PN0 (y)
N
X
L0j (y) f (xj ),
j=0
j = 0, . . . , N.
155
x.
= 0,
0 x0 + 1 x1 + + N xn = 1,
..
..
.
.
0 xk0 + 1 xk1 + + N xkn = ky k1 .
Para k = N , el sistema proporciona la misma regla de derivacion numerica que el metodo interpolatorio. Por u
ltimo, el metodo del desarrollo de Taylor requiere para este caso las etapas
siguientes:
(i) Desarrollar por Taylor, en un punto conveniente, los valores f (xj ) que aparecen en (5.7).
(ii) Desarrollar por Taylor en el mismo punto la derivada f 0 (y).
(iii) Igualar ambos desarrollos y deducir los coeficientes y el grado de exactitud, resolviendo el
sistema lineal resultante.
Ejemplo 3.7. Con N = 1, el polinomio P1 tiene la forma
P1 (x) = f (x0 )L0 (x) + f (x1 )L1 (x);
L0 (x) =
156
(x x1 )
,
x0 x1
L1 (x) =
(x x0 )
.
x1 x0
Figura 5.4: Comparacion de errores con evaluaciones del integrando para f (x) =
La aplicacion del metodo interpolatorio proporciona los coeficientes
0 =
1
,
x0 x1
1 =
1
.
x1 x0
y la formula
D2 (f ) =
f (x1 ) f (x0 )
.
x1 x0
f (y + h) f (y)
,
h
h = x1 x0 .
f (y + h) f (y h)
,
2h
x.
5.5.1.
Aproximaci
on a derivadas de orden superior
Los mismos procedimientos pueden servir para aproximar evaluaciones de derivadas de orden
superior. Se muestran los siguientes ejemplos como ilustracion [6, 5, 7]
Ejemplo 3.8. La aplicacion del metodo directo para aproximar la derivada segunda en un
punto y mediante una formula del tipo
f 00 (y) 0 f (x0 ) + 1 f (x1 ) + 2 f (x2 ),
con x0 = y h, x1 = y, x2 = y + h, h > 0, da lugar al sistema
0 + 1 + 2 = 0,
0 (y h) + 1 y + 2 (y + h) = 0,
0 (y h)2 + 1 y 2 + 2 (y + h)2 = 2;
de donde se obtienen los coeficientes 0 = 2 = 1/h2 , 1 = 2/h2 y un grado de exactitud
M 2. Los coeficientes obtenidos satisfacen la condicion de exactitud para f (x) = x3 :
0 (y h)3 + 1 y 3 + 2 (y + h)3 = 6y,
pero no la condicion para f (x) = x4 , por lo que M = 3.
Ejemplo 3.9 [6]. La aproximacion a la derivada tercera con una formula de cinco puntos
f 000 (y) 0 f (x0 ) + 1 f (x1 ) + 2 f (x2 ) + 3 f (x3 ) + 4 f (x4 ),
con x0 = y 2h, x1 = y h, x2 = y, x3 = y + h, x4 = y + 2h, h > 0, proporciona por cualquiera de
los tres metodos, los coeficientes 0 = 1/(2h3 ), 1 = 1/h3 , 2 = 0, 3 = 1/h3 , 4 = 1/(2h3 ) y
un grado de precision M = 4.
Las tablas 4.1 y 4.2 de [9] muestran los coeficientes de otras formulas de aproximacion.
5.5.2.
Error en la derivaci
on num
erica
f M +1) (t)K(t)dt,
donde K es la funci
on llamada n
ucleo de Peano, cuya expresi
on es
K(t) =
1
E(gt ),
M!
158
siendo gt la funci
on
gt (x) =
(x t)M
0
si
si
xt
x<t
Si, adem
as, el n
ucleo K no cambia de signo en [a, b], entonces para cada f C M +1 ([a, b]) existe
[a, b] tal que
E(f ) =
f M +1) ()
E(xM +1 ).
(M + 1)!
Demostraci
on. Es similar a la del teorema 5.3.1. Usando tambien la formula de Taylor con
resto integral
Z b
1
(x a)M M )
f M +1) (t)gt (x)dt,
f (a) +
f (x) = f (a) + +
M!
M! a
aplicando E a ambos lados, usando la hipotesis de que el grado de precision es al menos M y la
derivacion bajo signo integral, se tiene
Z b
Z b
1
1
M +1)
M +1)
E(f ) =
E
f
(t)gt (x)dt =
f
(t)E(gt )(x)dt .
M!
M!
a
a
La segunda parte del teorema se demuestra igual que el Corolario 5.3.2.
Ejemplo 3.10. Buscamos expresion del error para la formula centrada de aproximacion a la
derivada primera
f 0 (y)
f (y + h) f (y h)
,
2h
h > 0.
El n
ucleo de Peano es
(
K(t) =
(ty+h)
h
2
(tyh)
2
si
yhty
si
y ty+h
Como el n
ucleo no cambia de signo en [y h, y + h], existe [y h, y + h] tal que si f
C 3 ([y h, y + h]),
E(f ) =
(2h)2 000
h2
f () = f 000 ().
24
6
5.5.3.
El car
acter de mal puesto de la derivaci
on num
erica
Las formulas de derivacion numerica son especialmente sensibles a los errores de redondeo
cometidos por una deficiente evaluacion de la funcion. Este comportamiento es consecuencia de
que la evaluacion de la derivada de una funcion f en un punto es lo que se llama un problema
mal puesto. Esto significa que cambios peque
nos en los valores de f pueden generar grandes
cambios en el valor de la derivada en un punto. El ejemplo mostrado en [6] (vease tambien [7])
159
es suficientemente ilustrativo. Para una funcion f que tome valores cuyo valor absoluto sea del
orden de uno, la funcion g(x) = f (x) + 1012 sin(1018 x) toma valores parecidos a los de f , pero
no as la derivada, pues, por ejemplo, g 0 (0) = f 0 (0) + 106 .
La influencia de esto en la aproximacion numerica puede ilustrarse con la formula centrada
para aproximar f 0 (0) [6]
f 0 (0)
f (h) f (h)
,
2h
h > 0,
para cierto . La formula nos dice que para que el error sea peque
no, ha de serlo tambien h.
Pero, entonces, la formula puede magnificar los errores cometidos en f (h) al multiplicarlos por
un factor h1 grande.
h
100
101
102
103
104
105
106
f (h)
2,718282
1,105171
1,010050
1,001001
1,000100
1,000010
1,000001
(h)
|f 0 (0) f (h)f
|
2h
0,175201
0,001668
0,000017
0,000017
0,000166
0,001358
0,016523
5.6.
Ejercicios
Ejercicio 1. Halla el n
ucleo de Peano del error para las reglas de los trapecios y la de Simpson.
Obten, a partir de el, estimaciones del error de cuadratura correspondiente.
Ejercicio 2. Elabora un pseudocodigo para implementar las reglas del rectangulo, del punto
medio, la regla de los trapecios y la de Simpson compuestas. El pseudocodigo debe tener como
datos de entrada el diametro de la particion (que, por tanto, se supone uniforme) y los extremos
del intervalo de integracion.
Ejercicio 3. (i) Se va a calcular la integral de exp(-x2 ) entre 0 y 1 por la regla del punto medio
compuesta usando una particion uniforme. Se desean errores menores que 5 107 . Cu
antos
subintervalos hay que usar?
(ii) Determina el valor de h requerido para aproximar
Z
e2x sin(3x)dx,
160
con un error de 104 usando: (a) la regla del punto medio compuesta, (b) la regla de los trapecios
compuesta, (c) la regla de Simpson compuesta.
el n
Ejercicio 4. Sea K el n
ucleo de Peano de la regla del punto medio (simple) y sea K
ucleo
= K 0 .
que se obtiene para la misma regla tomando N = 0. Integrando por partes prueba que K
Cambia de signo K?
Ejercicio 5. Considerese una formula de cuadratura de grado de precision lo mayor posible de
la forma
Z b
ba
f (x)dx f (a) + f (a +
) + f (b),
4
a
para un intervalo cualquiera [a, b]. Calcula el n
ucleo de Peano correspondiente, y encuentra una
estimacion para el error tanto cuando se aplica la regla de cuadratura en forma simple como
cuando se aplica en forma compuesta.
Ejercicio 6. Considera una regla de cuadratura
Z b
f (x)dx 0 f (x0 ) + 1 f (x1 ) + + N f (xN )
a
que integra exactamente funciones constantes. Demuestra que el grado de precision de la formula
es M si y solo si para cualquier n
umero real la regla no da error al cuadrar las funciones
1, (x ), (x )2 , . . . , (x )M y da error no nulo al cuadrar (x )M +1 . Cual sera una buena
eleccion de para hacer menos cuentas al calcular el grado de precision?
Ejercicio 7. Regla de los trapecios modificada. Tambien pueden usarse reglas de cuadratura
que requieran conocer derivadas de f ; por ejemplo
Z b
f (x)dx 0 f (a) + 0 f 0 (a) + 1 f (b) + 1 f 0 (b).
a
son interpolatorias y que los nodos, todos distintos, satisfacen las relaciones
xj = a +
ba
(yj c),
dc
161
0 j N.
Demuestra que
j =
ba
j ,
dc
0 j N.
tenga grado de precision maximo. Generaliza dicha formula de tal manera que aproxime la integral
en un intervalo dado [a, b] cualquiera.
Ejercicio 10. Cuadratura gaussiana. Aparte de los pesos, en una regla de cuadratura tambien
se pueden elegir los nodos. Determina 0 , 1 , x0 , x1 de manera que la formula
Z 1
f (x)dx 0 f (x0 ) + 1 f (x1 ),
1
que sea exacta para polinomios cuadraticos. Aplica la formula para estimar
Z 0,1
ln(2 + x)dx.
0
Ejercicio 12. Estudie la convergencia para las formulas progresiva y regresiva de aproximaci
on
a la derivada de una funcion en un punto.
Ejercicio 13. Enuncie y demuestre la version del teorema de Peano para una formula de aproximacion a f 00 (y) del tipo (5.7). Estudie con ella la convergencia de la aproximacion
f 00 (0)
h > 0.
Ejercicio 14. Encuentre una formula que aproxime f 00 (0) usando una combinacion de los valores
f (0), f 0 (0), f (h), f 0 (h), utilizando el metodo interpolatorio y el metodo directo. Calcule el n
ucleo
de Peano del error y estudie la convergencia de la aproximacion cuando h 0.
162
Referencias
163
Tema 6
Resoluci
on num
erica de ecuaciones
diferenciales
1. Introduccion.
2. Metodos Runge Kutta para PVI en EDOs.
3. Metodos numericos para problemas de contorno.
4. Resolucion numerica de EDPs de evolucion.
5. Ejercicios.
6.1.
Introducci
on
A partir del nacimiento del Calculo Infinitesimal en el siglo XVII, surgieron ciertos problemas,
en Geometra y Mecanica, que dieron lugar a ecuaciones diferenciales de primer orden. A ello le
siguio un estudio sistematico de las ecuaciones diferenciales en el siglo XVIII y particularmente en
el XIX, como consecuencia de la aplicacion de las leyes de Newton de la Mecanica, para determinar
las trayectorias de cuerpos celestes. El hecho de que solo un n
umero limitado de ecuaciones
pueda resolverse por cuadratura y que metodos como el de desarrollo en serie (que parte de una
idea de Newton) no fueran del todo satisfactorios por dar lugar a representaciones infinitas de
las soluciones, todo ello explica el interes en la construccion de metodos de aproximaci
on. El
primero, usado por Euler en el siglo XVIII, para encontrar una aproximacion poligonal a la curva
solucion, dio lugar a diferentes mejoras en el siglo XIX, cuyo estmulo proceda de los metodos
de cuadratura. As, el metodo de diferencias finitas que lleva a la formula de Gregory tiene su
correspondiente aplicacion a ecuaciones diferenciales en los trabajos de Adams. Las reglas del
trapecio, Simpson y las gaussianas fueron la inspiracion para los trabajos de Runge, Heun y
Kutta. Estos u
ltimos metodos son hoy conocidos como metodos de un paso, mientras que la
tecnica propuesta por Adams es un metodo multipaso.
El uso de estos esquemas es hoy esencial, en parte porque cada vez mas disciplinas requieren
el uso de ecuaciones diferenciales para sus modelos matematicos y en parte porque estas tecnicas
numericas, basadas en la discretizacion de la ecuacion diferencial, son particularmente adaptables
164
al ordenador. Esta
es la razon por la que estos metodos, conocidos desde hace mas de un siglo,
solo hayan tenido un pleno desarrollo a partir de los a
nos cincuenta.
Este tema intenta ser una introduccion a la resolucion numerica de ecuaciones diferenciales.
Queda dividido en tres partes, en funcion del problema que se trata. las dos primeras abordan,
respectivamente, la aproximacion a problemas de valores iniciales y a problemas de contorno,
referidos a ecuaciones diferenciales ordinarias (EDOs, de ahora en adelante). la tercera parte describe metodos elementales de integracion numerica de ecuaciones en derivadas parciales (EDPs,
de ahora en adelante) de evolucion. Dado el caracter introductorio, y por tanto limitado, de la
leccion, en esta parte se consideran las ecuaciones del calor y de ondas como modelos a estudio
e ilustracion de los procedimientos numericos.
Como hemos indicado, la necesidad de metodos numericos para estudiar este tipo de problemas
procede del hecho de que, en general, una ecuacion diferencial no se puede integrar; no es posible
encontrar expresiones analticas de las soluciones. Tambien ocurre a veces que, aunque se tengan
expresiones teoricas, estas no sean manejables en la practica (como series infinitas, por ejemplo).
Como ilustracion de todo esto, se puede mostrar el modelo del circuito no lineal siguiente, en el
que la resistencia depende de la intensidad, y tiene por ecuacion
(1 +
2
1
1 dE
di
i(t) =
i(t)) +
R0
dt R0 C
R0 dt
C Condensador
@
@
i(t) Intensidad
R
E(t) Batera
R(t) = R0 + i(t)
E(t)
i(t)
165
(1 +
2
di
1
1 dE
i(t)) +
i(t) =
R0
dt R0 C
R0 dt
M
etodo num
erico
Aproximaciones i1 , i2 , . . . , iN a i(t1 ), i(t2 ), . . . , i(tN ) en instantes de tiempo t1 , t2 , . . . , tN
ik i(tk ),
k = 1, 2, . . . , N
6.2.
M
etodos Runge Kutta para PVI en EDOs
Esta seccion se dedica a intentar describir algunos de los metodos utilizados para tal fin en
el caso de un problema de valores iniciales para un sistemas de EDOs y, como toda la lecci
on,
tiene un caracter introductorio. Su estructura es la siguiente. Mencionaremos primero algunas
ideas, a tener en cuenta en el resto de la seccion, referentes al problema continuo de valores
iniciales a aproximar, as como a los metodos de aproximacion. Los primeros ejemplos de estos nos
serviran como introduccion para los llamados metodos generales de un paso. Luego, trataremos
166
6.2.1.
Conceptos generales
t0 t T,
(6.1)
(6.2)
aproximaciones. Este
es el significado de la siguiente propiedad: sean (t) C([a, b]), 0 RN y
consideremos el problema perturbado asociado a (6.1):
d
z = f (t, z) + (t),
dt
z(t0 ) = y0 + 0 .
t0 t T,
Se dice que el problema (6.1) esta bien puesto si se verifican las condiciones:
167
(6.3)
n = 0, . . . , N,
h=
tN t0
.
N
Nuestra descripcion, por sencillez, utilizara longitud de paso constante h. El uso de longitudes de
paso variables tiene que ver con la dificultad de aproximacion numerica a la solucion de (6.1),
de manera que se suele adaptar la red en funcion de la variacion de esta, si es conocida. La regla
general dice que suele ser mas apropiado utilizar longitudes hn mas peque
nas cuando la soluci
on
cambia rapidamente (por ejemplo, si oscila mucho) y longitudes mayores cuando esa variaci
on es
mas lenta [8].
168
n = 0, 1, . . . , N 1,
(6.4)
genera el llamado metodo de Euler explcito. La formula (6.4) procede de aproximar la derivada
en tn por un cociente incremental progresivo entre tn+1 y tn ,
dy
y(tn+1 ) y(tn )
y(tn+1 ) y(tn )
=
|
,
dt tn
tn+1 tn
h
lo que geometricamente significa aproximar la recta tangente a y en tn por la recta secante que
pasa por (tn , yn ), (tn+1 , yn+1 ) (figura 6.2).
n = 0, 1, . . . , N 1.
(6.5)
A diferencia de (6.4), la ecuacion (6.5) requiere resolver en cada paso un sistema no lineal para
obtener yn+1 con alg
un procedimiento iterativo, es decir, es lo que se llama un metodo implcito.
169
.
dt t=tn+1
h
Ejemplo 8.3. Metodo de Heun. En los ejemplos anteriores, el avance de la solucion se produce
seg
un una pendiente calculada en cada paso. Una alternativa viene dada mediante la sustituci
on
de este valor por alguna pendiente media de y en el intervalo [tn , tn+1 ]. Esta idea generar
a la
llamada clase de los metodos Runge-Kutta, que luego trataremos. Un elemento de esta familia
es el llamado metodo de Heun. Aqu la solucion avanza en cada paso seg
un el promedio de tres
pendientes: la dada en el punto (tn , yn ) y las dadas en los puntos (tn + h/3, yn + h3 f (tn , yn )) y
h
(tn + 2h/3, yn + 2h
ormula de
3 f (tn , yn+1/3 )), donde yn+1/3 = yn + 3 f (tn , yn ). De esta manera, la f
avance queda
yn+1 = yn +
h
(k1 + 3k3 ),
4
n = 0, 1, . . . , N 1,
(6.6)
donde
k1 = f (tn , yn ),
h
k1 ),
3
2h
= f (tn + 2h/3, yn +
k2 ).
3
k2 = f (tn + h/3, yn +
k3
A pesar de ser mas larga, la formula (6.6) sigue siendo explcita y de un paso. La diferencia con las
anteriores esta en el n
umero de pendientes aproximativas que se utilizan para definir el metodo,
tambien llamado n
umero de etapas, que en este caso es tres.
Ejemplo 8.4. Este procedimiento puede tambien generar metodos implcitos de varias etapas.
Por ejemplo,
yn+1 = yn +
h
(k1 + k2 ),
2
n = 0, 1, . . . , N 1,
donde
k1 = f (tn , yn ),
k2 = f (tn + h, yn +
h
h
k1 + k2 ).
2
2
n 0,
sean peque
nos para todos los tn . Sin embargo, los esquemas no intentan controlar estas cantidades. Mas bien incorporan informacion de la ecuacion diferencial en los k pasos anteriores
(tn , yn ), . . . , (tnk+1 , ynk+1 ). En particular, los metodos de un paso tratan de aproximar de forma optima la solucion local un (t) o solucion del problema
u0 (t) = f (t, u),
u(tn ) = yn .
De esta manera, el error global puede verse como la suma de los diferentes errores locales (figura
6.3)
eln+1 = u(tn+1 ) yn+1 ,
n 0,
El orden de un metodo mide la magnitud del error local que el metodo comete en un paso, de
manera que se dice que un metodo es de orden p si los errores locales cumplen
||u(tn + h) yn+1 || Chp+1 ,
con C una constante. Si p 1, entonces cuanto menor es h (es decir, cuanto mas densa es la red)
menor es el error local en tn+1 . Mas adelante nos ocuparemos de estas cuestiones.
171
6.2.2.
M
etodos de un paso
Nos centramos ahora en los metodos de un paso. La estrategia sera presentar primero los
llamados metodos de Taylor, que son los metodos de un paso conceptualmente mas sencillos
(estan basados en desarrollos de Taylor de la solucion) y que nos serviran para introducir todos
los elementos para la descripcion general que realizaremos despues.
M
etodos de Taylor
Por sencillez, describiremos los metodos para una ecuacion escalar, d = 1 en (6.1). La extensi
on
a un sistema de EDOs es directa, aunque los calculos son mas incomodos.
La idea basica de los metodos de Taylor es la siguiente: si la solucion y(t) es suficientemente
regular, podemos desarrollar en serie de Taylor
y(t + h) = y(t) + hy 0 (t) + +
hq q)
y (t) +
q!
(6.7)
hp p)
y .
p! n
(6.8)
j)
donde yn denota la derivada j-esima evaluada en (tn , yn ). La formula (6.8) puede escribirse de
manera mas formal
yn+1 = yn + h(tn , yn , h),
(6.9)
con
(t, y, h) = y 0 +
hp1 p)
h 00
y + +
y .
2
p!
Estos metodos son relativamente efectivos si las derivadas totales del segundo miembro de (6.1)
no son difciles de evaluar (hay software de calculo simbolico disponible para ello: MATLAB,
MAPLE, MATHEMATICA, etc). Aunque poco utilizados, en comparacion con otros metodos
de un paso, la potencia actual de los ordenadores, unida a la efectividad del paquete simb
olico
utilizado, hace que los metodos de Taylor no sean totalmente descartables. En cualquier caso,
suele ser generalmente mas efectivo tratar de buscar metodos de un paso que no precisen de
tales evaluaciones explcitas de las derivadas. De esta idea surgiran los metodos Runge-Kutta,
que trataremos posteriormente.
Ejemplo 8.5. Ejemplo de un metodo de Taylor. Podemos usar el metodo de Taylor de orden
p = 3 para aproximar el problema de valores iniciales
d
ty
y(t) = f (t, y) =
,
dt
2
y(0) = 1.
172
t [0, 3],
y 0 (t) = f (t, y) =
y 00 (t) =
y 000 (t) =
la formula del metodo queda:
h3
h2
yn+1 = yn + hyn0 + yn00 + yn000
2 6
tn y n
h 2 2 tn + y n
h3 2 + tn yn
= yn + h
+
+
,
2
2
4
6
8
donde tn = nh,
n = 0, 1, 2, . . ..
Forma general de un m
etodo explcito de un paso
La formula (6.9) puede utilizarse para introducir metodos mas generales de un paso, explcitos,
cuya expresion dependera de la eleccion de la funcion incremento . Se imponen dos hip
otesis
sobre esta funcion:
(t, y, 0) = f (t, y),
y ||,
||(t, y, h) (t, y , h)|| L||y
(6.10)
(6.11)
n = 0, 1, . . . , N,
173
se llaman errores globales. La magnitud de la diferencia entre la solucion exacta en los nodos
{y(tn )}N
erica {yn }N
n=0 y la num
n=0 se define como
max ||y(tn ) yn ||,
0nN
max ||y(tn ) yn || = 0,
h0+ 0nN
cuando y0 (h) y(t0 ), h 0. Se dice que la convergencia es de orden p si, supuesto que y0 (h)
y(t0 ) cuando h 0, p es el maximo entero tal que
max ||y(tn ) yn || = O(hp )
(h 0+),
0nN
t > tn ,
u(tn ) = yn ,
se llama solucion local de (6.1) en (tn , yn ), mientras que
eln+1 = u(tn+1 ) yn+1 = u(tn+1 ) yn h(tn , yn , h),
(6.12)
se llama error local del metodo en tn+1 . Hay que observar que el error local compara la soluci
on
numerica en tn+1 con la solucion del problema desde (tn , yn ). Se dice entonces que el metodo es
de orden p si los errores locales verifican
||eln+1 || = O(hp+1 ),
h 0+,
n = 0, 1, . . . , N,
para cualquier problema (6.1) con f C p verificando las hipotesis del teorema 6.2.1. En la
practica, el orden del error local se determina mediante convenientes desarrollos de Taylor de
(6.12). De esta manera, el orden p equivale a errores locales O(hp+1 ) o a errores globales O(hp ).
Error local de truncaci
on y consistencia
Dado que la definicion de la convergencia de un metodo involucra a la solucion teorica, por
regla general desconocida, debe buscarse un procedimiento alternativo e indirecto para el estudio
174
(6.13)
Es decir, el error local cometido en tn+1 , si se parte de la solucion teorica y(tn ) en tn . Observemos
que, si yn es una buena aproximacion a y(tn ), entonces
yn+1 = yn + h(tn , yn , h) y(tn ) + h(tn , y(tn ), h).
Como queremos que yn+1 sea una buena aproximacion a y(tn+1 ), necesariamente
y(tn+1 ) yn+1 y(tn+1 ) y(tn ) h(tn , y(tn ), h),
ha de ser peque
no. Esta
es la razon de la introduccion del error local de truncacion.
Se dice que el metodo (6.9) es consistente si para todo problema de valores iniciales (6.1), con
f verificando las hipotesis del teorema 6.2.1, se tiene
eltn+1
= 0.
lm
max
h0+ 0nN 1
h
Se dice que la consistencia es de orden p si ademas, cuando f C p ,
eltn+1
= O(hp ), h 0 + .
max
0nN 1
h
Ejemplo 8.8. Veamos que el metodo de Euler explcito (6.4) es consistente de orden uno.
Observemos que
eltn+1 = y(tn+1 ) y(tn ) hf (tn , y(tn ))
= y(tn+1 ) y(tn ) hy 0 (tn ).
Desarrollando por Taylor en el punto (tn , y(tn )) y suponiendo que y C 2 (de manera que f C 1 )
eltn+1 = (y(tn ) + hy 0 (tn ) +
h2 00
h2 00
y (n )) y(tn ) hy 0 (tn ) =
y (n ),
2
2
(6.14)
de manera que
eltn+1
= O(h),
max
0nN 1
h
h 0+,
175
Estabilidad
Para obtener unos errores globales peque
nos, no basta con buscar errores de truncaci
on
peque
nos, pues al hacer h 0, estos errores pueden disminuir, pero, al mismo tiempo, n aumenta, lo que puede generar una propagacion de errores incontrolada, al aumentar el n
umero de
operaciones. Para analizar esto se introduce el concepto de estabilidad.
N
Observemos que las sucesiones {yn }N
n=0 e {y(tn )}n=0 satisfacen, respectivamente, las ecuaciones en diferencias
(6.15)
(6.16)
con valores de arranque y0 y z0 = y0 + d0 , y donde las perturbaciones {dn }n vienen dadas por
los errores de truncacion local.
Se dice que el metodo (6.9) es estable si para todo problema (6.1) con f verificando las
hipotesis del teorema 6.2.1, las soluciones de (6.15), (6.16) verifican
max ||yn zn || S
0nN
N
X
||dj ||,
j=0
con S > 0 constante (llamada constante de estabilidad), independiente de h y que solo depende
de los datos de (6.1).
Esta propiedad representa el hecho de que peque
nas perturbaciones del esquema en diferencias
que define el metodo, medidas en la norma
||d|| =
N
X
d = {dj }N
j=0 ,
||dj ||,
j=0
(1 + hL)N ||d0 || +
n
X
j=1
176
N
X
||dj ||,
(6.17)
j=0
con S = eL(T t0 ) .
Busquemos ahora la relacion entre las condiciones (6.10), (6.11) con las dos propiedades que
acabamos de definir, consistencia y estabilidad. Por un lado, supuesto que verifica (6.10) y es
suficientemente regular, se tiene que
(t, y, h) = (t, y, 0) + hh (t, y, 0) + O(h2 )
= f (t, y) + hh (t, y, 0) + O(h2 ),
h0+.
Entonces
eltn+1 = y(tn+1 ) y(tn ) h(tn , y(tn ), h)
= y(tn+1 ) y(tn ) hf (tn , y(tn )) + O(h2 ),
h 0+,
de manera que el metodo es consistente (con orden al menos uno). Luego la condicion (6.10)
implica consistencia. Por otro lado, la relacion entre (6.11) y la estabilidad viene dada a traves
de un razonamiento similar al realizado en el ejemplo 8.8. Si satisface (6.11),
||yn zn || = ||yn1 zn1 + h((tn1 , yn1 , h) (tn1 , zn1 , h)) + dn ||
n1 zn1 || + ||dn ||
||yn1 zn1 || + hL||y
n1 zn1 || + ||dn ||.
= (1 + hL)||y
Luego, al igual que antes, se tiene
||yn zn || S
N
X
||dj ||,
j=0
177
= e
n=0
2
L(T t0 )
MNh = e
(T t0 )M h,
0nN
h 0,
6.2.3.
La clase de los m
etodos Runge-Kutta
Formulaci
on
Al tratar los metodos de Taylor, comentamos la conveniencia de cambiar la estrategia del
calculo de derivadas de la solucion para obtener metodos numericos de un paso, por la dificultad
de computacion de las mismas. Una alternativa viene dada por los metodos Runge-Kutta (RK).
Estos
proceden de la idea de utilizar la formulacion integral del problema,
Z tn +h
u(tn + h) = u(tn ) +
f (t, u(t))dt,
tn
y una regla de cuadratura para aproximar la integral. Por ejemplo, el metodo de Euler explcito
surge cuando se usa la regla del rectangulo basada en el punto tn . Cuando se desarrolla dicha
formula en tn+1 , se obtiene el metodo de Euler implcito (6.5). La regla de los trapecios proporciona el metodo implcito conocido como formula de los trapecios
yn+1 = yn +
h
(f (tn , yn ) + f (tn+1 , yn+1 )) .
2
i=1
s
X
(6.18)
j=1
Los valores intermedios u(tn + ci h) pueden aproximarse usando de nuevo la formulacion integral
Z tn +ci h
u(tn + ci h) = u(tn ) +
f (t, u(t))dt,
tn
178
y una formula de cuadratura para la integral. Un aspecto muy importante es que estas aproximaciones intermedias Yi no requieren reglas de cuadraturas tan precisas, de modo que es posible
mantener un cierto orden alto en la formulacion (6.18) con aproximaciones Yi de orden menor.
Por ejemplo, utilizando la regla del punto medio en (6.18) y aproximando el valor intermedio
u(tn + h2 ) mediante el metodo de Euler explcito (es decir, aproximando la integral intermedia
por la regla del rectangulo basada en el punto tn ) se obtiene la llamada regla del punto medio
h
h
yn+1 = yn + hf tn + , yn + f (tn , yn ) .
2
2
Presentados estos ejemplos, a modo de introduccion, pasamos a la formulacion general de un
Yi = yn + h
s
X
i=1
i1
X
bi f (tn + ci h, Yi ),
aij f (tn + cj h, Yj ),
(6.19)
i = 1, . . . , s,
j=1
(con Y1 = yn ). Las cantidades Yi son las etapas internas del metodo y el esquema queda identificado por los coeficientes bi , ci , i = 1, . . . , s, aij , i = 1, . . . , s, j = 1, . . . , i 1. De este modo la
representacion de un metodo RK suele darse a traves del llamado tablero de Butcher
0
c2
..
.
0
a21
..
.
cs
as1
b1
0
bs
A
bT
con
0
a21
A=
...
as1
0
0
..
.
as2
0
0
..
.
as,s1
0
0
,
..
.
b = (b1 , . . . , bs )T ,
c = (0, c2 , . . . cs )T ,
aij = ci ,
i = 1, . . . , s.
j=1
1. Metodo de Euler.
0
0
1
2. Metodo de Euler modificado.
0
1/2
0
1/2
0
0
0
1
0
1
1/2
0
0
1/2
0
1/3
0
1/4
0
0
2/3
0
0
0
0
3/4
0
1/2
1
1/6
0
0
2
2/3
0
0
0
1/6
0
1/2
0
0
1/6
0
0
1/2
0
1/3
0
0
0
1
1/3
0
0
0
0
1/6
s
X
bi ki ,
i=1
ki = f tn + ci h, yn + h
i1
X
aij kj ,
i = 1, . . . , s,
j=1
s
X
bi f (tn + ci h, Yi ),
i=1
180
Yi ||
= ||y y + h
i1
X
aij f (tn + cj h, Yj ) f (tn + cj h, Yj ) ||
j=1
||y y || + h
i1
X
|aij |L||Yj Yj ||
j=1
i1
X
||y y || + h
|aij | L max ||Yj Yj ||
1ji1
j=1
1js
Yj ||
1
1 h||A|| L
||y y ||.
Entonces
s
X
bi (f (tn + ci h, Yi ) f (tn + ci h, Yi ))
i=1
s
X
|bi |L||Yj Yj ||
i=1
s
X
|bi |L
i=1
1
1 h||A|| L
!
||y y ||.
De manera que si h se elige con h |A|| L < 1, entonces para h [0, h ] el metodo es estable. En
realidad, para el caso de metodos RK explcitos, se tiene que todos ellos son estables. La raz
on
es la siguiente. Fijemonos en que
||Y1 Y1 || = ||y y ||
||Y2 Y2 || ||y y || + h|a21 |L||y y || (1 + h||A|| L)||y y || eh||A|| L ||y y ||
||Y3 Y3 || ||y y || + h (|a31 |L||y y || + |a32 |L||Y2 Y2 ||)
||y y || + h||A|| L (||y y || + ||Y2 Y2 ||)
(1 + h||A|| L)||y y || + h||A|| L(1 + h||A|| L)||y y ||
= (1 + h||A|| L)2 ||y y || e2h||A|| L ||y y ||.
As, suponiendo que
||Yj Yj || e(j1)h||A|| L ||y y ||,
entonces
||Yi Yi || ||y y || + h
i1
X
|aij |L||Yj Yj ||
j=1
181
j = 1, . . . , i 1,
||y y || + h
i1
X
j=1
||y y || + h
i1
X
j=1
Luego
s
X
bi (f (tn + ci h, Yi ) f (tn + ci h, Yi )) ||
i=1
s
X
|bi |L||Yj
Yj ||
i=1
s
X
!
|bi | Le(s1)h||A|| L ||y y ||
i=1
y ||,
= L||y
donde
=
L
s
X
!
|bi | Le(s1)h||A|| L .
i=1
i=1
bi = 1.
i=1
Introducci
on al estudio del orden de un m
etodo RK
Entonces, seg
un lo comentado previamente, un metodo RK explcito, de s etapas, es convergente si
s
X
bi = 1.
i=1
Otro aspecto a considerar es el orden de convergencia. Una medida del coste computacional de
un metodo RK viene dada por las evaluaciones de la funcion f de (6.1), que se supone aporta el
mayor gasto operativo. Puesto que un metodo RK explcito de s etapas requiere, para avanzar un
paso, s evaluaciones de f , es importante obtener un orden de convergencia alto con un n
umero
bajo de etapas.
La manera mas natural de construir formulas RK de un orden dado p consiste en imponer
directamente, sobre los coeficientes del tablero de Butcher, las condiciones derivadas de exigir
que
||y(tn + h) yn+1 || = O(hp+1 ),
182
h 0.
Esto equivale a requerir coincidencia, hasta terminos de orden p, de los desarrollos de Taylor
en h de y(tn + h) e yn+1 (h), con f suficientemente derivable, y suponiendo que yn+1 es el valor
obtenido por el metodo partiendo de la solucion exacta en tn , de modo que y(tn + h) yn+1 (h)
es el error de truncacion local en tn+1 . Debemos buscar entonces consistencia de orden p.
Como ilustracion del procedimiento consideremos que (6.1) es una ecuacion escalar (d = 1). Si
deseamos determinar metodos de una etapa, los desarrollos solo se pueden hacer coincidir hasta
terminos de primer orden si b1 = 1 pues
y(tn + h) = y(tn ) + hf (tn , y(tn )) + O(h2 ),
yn+1 (h) = y(tn ) + hb1 f (tn , y(tn )),
de modo que solo hay un metodo explcito de orden uno, el de Euler (6.4). Para metodos de dos
etapas
yn+1 = yn + h (b1 f (tn , Y1 ) + b2 f (tn + c2 h, Y2 )) ,
Y1 = yn
Y2 = yn + ha21 f (tn , Y1 ),
el desarrollo del metodo sera
yn+1 = y + h(b1 + b2 )f + h2 b2 a21 F + O(h3 ),
con y = y(tn ), f = f (tn , y(tn )), F = (ft + f fy )|(tn ,y(tn )) , mientras que el de la solucion es
y(tn + h) = y + hf +
h2
F + O(h3 ).
2
0
a21
a31
b1
0
a32
b2
0
b3
183
10
s Etapas
11
12 s 17
13 s 17
Tabla 6.1: N
umero de etapas necesario para algunos ordenes.
M
etodos RK generales
Los metodos RK no tienen por que ser explcitos. la formulacion general de un metodo RK
de s etapas para aproximar (6.1) tiene la forma
yn+1 = yn + h
Yi = yn + h
s
X
i=1
s
X
bi f (xn + ci h, Yi ),
aij f (xn + ci h, Yj ),
j=1
184
(6.20)
i = 1, . . . , s,
b2
A
,
=
bT
bs
En general, las formulas (6.20) corresponden a las de un metodo implcito porque el conjunto de
formulas que definen las etapas internas forman un sistema no lineal de s d ecuaciones en las
s d componentes de Yi , i = 1, . . . , s. Hay dos casos particulares destacables. El primero es el
de los metodos semiimplcitos. Aqu, A es triangular inferior con aii 6= 0 para alg
un ndice i. De
esta manera, el sistema de ecuaciones no lineales se desacopla en s sistemas, uno para cada etapa
interna Yi . El segundo caso particular es el de los metodos explcitos, ya comentados, en los que
A es estrictamente triangular inferior (con aii = 0 para todo i) y las etapas internas se obtienen
explcitamente.
Cuando el metodo no es explcito, se genera un problema a
nadido, referente a la existencia
N
de solucion numerica {yn }n=0 o sucesion que satisface la formula del esquema. Por otro lado,
el estudio del orden de los metodos es similar al explicado para los explcitos. En general, los
metodos RK implcitos consiguen, con el mismo n
umero de etapas que un explcito, orden mayor.
Ejemplo 8.10. Tableros de Butcher de varios metodos RK implcitos.
1. Metodo de Euler implcito.
1
1
1
2. Metodo de punto medio.
1/2
1/2
1
1
2
1
2
3
6
3
6
1
4
1
4
1
4
+
1
2
3
6
1
4
3
6
1
2
Implementaci
on de m
etodos de un paso
Vamos a ilustrar, a traves de un ejemplo, algunos aspectos de programacion de metodos
RK. Nos centraremos, por sencillez y dado el caracter preliminar de la leccion, en metodos
RK explcitos. Como temas que van mas alla de lo que aqu comentaremos se encuentran la
185
implementacion de metodos implcitos y el uso de paso variable [2, 4]. Son muy importantes, pero
quedan fuera del proposito introductorio de la leccion.
Consideremos, como metodo ilustrativo, la regla de Heun de tercer orden, con tablero
0
0
1/3
1/3
2/3
2/3
1/4
3/4
1
,
2x(t)
x(0) = 1,
1 + 3e8t
, v(t) = 3e8t , que se ha aproximado desde t = 0 hasta
8
t = 1 con h = 1/4, 1/8, 1/16, 1/32 y cuyos errores finales en t = 1 se muestran a la derecha de la
tabla 6.2.
h
0.5
0.25
0.125
0.0625
0.03125
0.015625
error
1.1971E-04
1.4273E-05
1.7268E-06
2.1186E-07
2.6222E-08
3.2612E-09
h
0.1
5E-02
2.5E-02
1.25E-02
6.25E-03
error
2.8635E-04
2.9445E-05
3.1707E-06
3.6621E-07
4.3984E-08
Tabla 6.2: Resultados de la regla de Heun para los dos problemas planteados.
Una forma de deducir el orden de un metodo a partir de una tabla de errores es la siguiente. Si
el metodo tiene orden p y E(h) denota el error obtenido con diametro h, entonces E(h) = O(hp ).
De manera que si h se divide por, digamos, dos, el error consiguiente E(h/2) debera dividirse, con
respecto a E(h), por 2p , aproximadamente. De esta manera, Se puede observar el orden tres del
metodo, si bien para el sistema es un poco mas alto. Las dos figuras 6.4, 6.5 muestran, en escala
logartmica, el error frente a la longitud de paso h para cada tabla. Las estrellas representan los
187
6.3.
6.3.1.
M
etodos num
ericos para problemas de contorno
Introducci
on
El segundo tipo de problemas cuya resolucion mumerica trataremos tiene la forma siguiente
y 00 = f (t, y, y 0 ),
y(a) = ,
a t b,
(6.21)
y(b) = .
a t b,
(6.22)
y(b) = .
6.3.2.
El m
etodo de tiro
El primer metodo numerico para (6.21) que consideraremos es el llamado metodo de tiro
simple. Esta tecnica involucra aproximaciones numericas a problemas de valores iniciales [3]. La
idea es conjeturar el valor inicial y 0 (a) para resolver, o aproximar, el problema de valor inicial
relacionado. Ello permite obtener una aproximacion a y(b). Si esta no coincide con el valor de
contorno , se puede corregir el valor de y 0 (a) conjeturado y comenzar de nuevo.
Caso lineal
Consideraremos primero el problema (6.22), del cual supondremos que admite solucion u
nica.
0
La conjetura y (a) = s1 (angulo de tiro) lleva al PVI
y 00 (t) + p(t)y 0 (t) + q(t)y(t) = r(t),
y(a) = ,
a t b,
y (a) = s1 ,
Denotemos por y (1) (t) a su solucion. En general, sera y (1) (b) 6= , por lo que podemos ensayar
un segundo angulo de tiro y 0 (a) = s2 e integrar el correspondiente PVI. Si denotamos por y (2) (t)
a su solucion, entonces la combinacion lineal
y(t) = y (1) (t) + (1 )y (2) (t),
con un parametro libre, es solucion de la ecuacion diferencial lineal y satisface la primera
condicion de contorno. Ahora, solo queda ajustar el valor de para que y(b) = , es decir
= y (1) (b) + (1 )y (2) (b) =
y (2) (b)
,
y (1) (b) y (2) (b)
siempre que y (1) (b) 6= y (2) (b). En la practica, los PVI son resueltos numericamente para calcu(1) (2)
lar aproximaciones yb , yb a y (1) (b), y (2) (b) respectivamente. Estos valores son utilizados para
calcular una aproximacion a
(2)
y( b)
(1)
(2)
yb yb
189
y 0 (a) = s1 + (1 )s2 ,
a t b,
(6.23)
y (a) = s.
denotamos por ys (t) = y(t, s) a su solucion. La idea es seleccionar el angulo de tiro s para que
y(b, s ) = o, lo que es lo mismo, resolver la ecuacion
(s) = y(b, s) = 0,
(6.24)
para encontrar s . El procedimiento requiere entonces, para cada s, resolver numericamente (6.23)
para obtener una aproximacion a y(b, s) y aplicar un metodo iterativo para resolver (6.24) con
dicha aproximacion, en lugar de y(b, s). Estos metodos iterativos (tema 4) exigen, por lo general,
evaluar en una sucesion de valores de s, con sus correspondientes aproximaciones a y(b, s). En
el caso de que f (t, u, v) sea lineal en u y v, entonces es lineal en s y bastan dos evaluaciones de
para resolver (6.24), que es lo que hicimos en el apartado anterior, al resolver el caso lineal. En
otro caso, el metodo requerira la aproximacion numerica de una sucesion de problemas de valor
inicial, cuyo n
umero dependera del n
umero de iteraciones del metodo utilizado para (6.24), por
lo que la eficiencia de este es importante. Por este motivo, un metodo clasico utilizado en este
contexto es el metodo de Newton, que adquiere aqu la forma
sn+1 = sn
(sn )
,
0 (sn )
n = 0, 1, . . . ,
con s0 dado. Para aproximar 0 (s), se deriva parcialmente (6.23) con respecto a s,
ys00
f t
f ys
f y 0
(t) =
+
+ 0 s,
s
t s ys s
ys s
0
ys
ys
(a) = 0,
(a) = 1,
s
s
que, si v = ys /s, se puede escribir
v 00 = fys (t, ys , ys0 )v + fys0 (t, ys , ys0 )v 0 ,
v(a) = 0,
(6.25)
v (a) = 1.
Observese que entonces 0 (s) = v(b). La aproximacion numerica al problema (6.25) proporciona
entonces una aproximacion a v(b), que se toma como valor de 0 (s) para el metodo de Newton.
190
M
etodo de tiro m
ultiple
Los principales inconvenientes del metodo de tiro m
ultiple tienen que ver con la longitud
del intervalo de integracion de los problemas de valor inicial involucrados, pues, como sabemos,
un intervalo de peque
no tama
no es mas conveniente, tanto para la existencia de solucion del
problemna teorico como para la estabilidad de la correspondiente aproximacion numerica.
Un procedimiento alternativo, denominado metodo de tiro m
ultiple, consiste en dividir el
intervalo [a, b] seg
un una particion
a = t0 < t1 < < tN = b
y aplicar el metodo de tiro simple en cada subintervalo [tj , tj+1 ], j = 0, . . . , N 1. Los valores
iniciales de cada PVI se determinan de modo que la solucion asociada tenga una regularidad C 1
en toda la particion, para obtener as una solucion de (6.21).
Una manera de implementar esta idea podra ser la siguiente: denotamos por yj (t, j , sj ), tj
t tj+1 a la solucion de
y 00 = f (t, y, y 0 ),
y(tj ) = j ,
tj t tj+1 ,
y (tj ) = sj ,
(6.26)
j = 0, . . . , N 1.
Si 0 = , la funcion
y(t, s0 , 1 , s1 , . . . , N 1 , sN 1 ) = yj (t, j , sj ),
si
tj t tj+1 ,
yN 1 (b, N 1 , sN 1 ) = ,
sistema que debe resolverse con un metodo iterativo. Una de las ventajas computacionales de este
procedimiento es que la integracion numerica de los problemas (6.26) puede hacerse en paralelo.
6.3.3.
M
etodo de diferencias finitas
Una alternativa al metodo de tiro viene dado por el llamado metodo de diferencias finitas.
Los pasos generales de este procedimiento son los siguientes. Primero, se discretiza el problema,
dividiendose el intervalo [a, b] en una red discreta de nodos a = t0 < t1 < < tN = b y
buscando aproximaciones yj a los valores de la solucion de (6.21) en el nodo tj , y(tj ). Tales
aproximaciones numericas se obtienen como solucion de una ecuacion en diferencias. A su vez,
esta surge de sustituir, en la ecuacion de (6.21) evaluada en tj , las derivadas por apropiadas
formulas de aproximacion (vease tema 5).
Caso lineal
El ejemplo mas utilizado para ilustrar este procedimiento esta asociado al problema lineal
(6.22) como modelo. La discretizacion considera, por sencillez, una red uniforme de nodos tj =
a + jh, j = 0, . . . , N , con h = (b a)/N , y las formulas en diferencias centradas para aproximar
las derivadas,
y 0 (tj ) =
y 00 (tj ) =
y(tj+1 ) y(tj1 )
+ O(h2 ),
2h
y(tj+1 ) 2y(tj ) + y(tj1 )
+ O(h2 ).
h2
191
(6.27)
(6.28)
con
h
h
pj , bj = 2 + h2 qj , cj = 1 + pj , j = 1, . . . , N 1,
2
2
2
2
dj = h rj , j = 2, . . . , N 2, d1 = h r1 a1 , dN = h2 rN 1 cN 1 .
aj = 1
b1 c1 0
0
0
0
a2 b2 c2
d1
0
0
0
d2
0 a3 b3
c3
0
0
..
..
.. . .
.
..
..
..
.
.
.
.
.
.
.
A= .
, d = . .
..
0 0 aN 4 bN 4 cN 4
0
0
0 0
dN 2
0
a
b
c
0
N
3
N
3
N
3
0 0
0
0
aN 2 bN 2 cN 2
dN 1
0 0
0
0
0
aN 1 bN 1
Algoritmo de Thomas
La resolucion numerica mas habitual de sistemas tridiagonales como el anterior viene dada por
una reformulacion compacta de la eliminacion gaussiana sin pivotaje para matrices tridiagonales,
conocida como algoritmo de Thomas (tema 2, ejercicio 6). La idea es reducir el sistema (6.28) a
otro de la forma
yj + j yj+1 = zj ,
j = 1, . . . , N 1,
(6.29)
para ciertos j , zj . Supuesto que se tiene (6.29) para j = 1, . . . k, se elimina yk de las ecuaciones
yk + k yk+1 = zk ,
ak+1 yk + bk+1 yk+1 + ck+1 yk+2 = dk+1 ,
obteniendose
yk+1 +
dk+1 ak+1 zk
ck+1
yk+2 =
.
bk+1 ak+1 k
bk+1 ak+1 k
cj
,
bj aj j1
zj =
dj aj zj1
, j = 1, . . . , N 1,
bj aj j1
192
(6.30)
= ,
= zj j yj+1 ,
j = N 1, . . . , 1.
(6.31)
Las formulas (6.30) estan bien definidas mientras el denominador no se anule. Debemos asegurar
ademas que la matriz A sea regular (para que el sistema tenga solucion) y que el algoritmo sea
estable, es decir, que los errores iniciales en yN no se amplifiquen a partir de la formula (6.31).
Puede comprobarse [5] que todo ello se consigue bajo la hipotesis de que
|b1 | > |c1 |,
es decir, que la matriz A sea estrictamente diagonalmente dominante. Ello impone restricciones
al valor del diametro h de la particion.
Convergencia
Es razonable pensar que las aproximaciones numericas a (6.22) as construidas se aproximen
cada vez mas a la solucion teorica en los nodos a medida que h se hace mas peque
no, pues la
ecuacion en diferencias procede de truncaciones en las formulas de las derivadas y cuyos errores
tienden a cero como h2 . En la seccion anterior, necesitabamos que estos errores no se propagaran de forma perjudicial para obtener errores globales aceptables. Es decir, requeramos que
las aproximaciones fueran estables, propiedad que, unida a la consistencia (errores de truncaci
on
peque
nos), nos llevaba a la convergencia (errores globales peque
nos).
La situacion en este caso es algo similar. Se define el error de truncacion j en el nodo
tj , j = 1, . . . , N 1, como el residuo obtenido al sustituir, en la ecuacion en diferencias (6.27), la
aproximacion yj por la solucion teorica y(tj ), esto es
j (h) :=
que tiene claramente que ver con la truncacion de las aproximaciones a las derivadas. Se dice que
el esquema en diferencias es consistente si j 0, j = 1, . . . , N 1 cuando h 0. Utilizando
las formulas de aproximacion a las derivadas, se tiene que para y(t) suficientemente regular,
j = O(h2 ), h 0, j = 1, . . . , N 1, y la consistencia de (6.27) es de orden dos.
Si definimos el error global en tj como ej = y(tj ) yj , j = 0, . . . , N , entonces ej satisface la
ecuacion en diferencias
193
max
1iN 1
|(Lh V )i |}},
max
1iN 1
|i |,
j = 1, . . . , N 1,
0 t /2,
y(/2) = 1,
error
2.7103E-03
6.6667E-04
1.6916E-04
4.2275E-05
1.0568E-05
6.4.
Resoluci
on num
erica de EDPs de evoluci
on
6.4.1.
Problemas parab
olicos
El modelo de ecuacion parabolica que manejaremos es la llamada ecuacion del calor, cuya utilizacion para estudiar (al menos, en una primera aproximacion) fenomenos de difusion es conocida.
Una buena discusion sobre la deduccion de la ecuacion y los distintos tipos de condiciones de
contorno puede consultarse en referencias como [1, 10, 6]. Por sencillez, nuestro problema modelo
incorporara una condicion de contorno de tipo Dirichlet
ut (x, t) = Duxx (x, t) + f (x, t),
u(x, 0) = u0 (x),
u(0, t) = g0 (t),
0 < x < L, 0 t T,
0 < x < L,
(6.32)
u(L, t) = gL (t),
donde
1. u = u(x, t) puede representar la tremperatura, en el instante t, en la posicion x de una
barra de longitud L, suficientemente fina como para considerarla unidimensional, de secci
on
transversal constante y de material termoconductor homogeneo. Tambien se supone que la
generacion de calor en el material de la barra es despreciable.
u(0, t) = g0 (t)
f
u(x, t)
L
u(L, t) = gL (t)
f
-
0 t T,
0 t T,
(6.33)
195
que
h = L/J, k = T /N,
(6.34)
xj = jh, j = 0, . . . , J,
tn = nk, n = 0, . . . , N.
El segundo paso consiste en aproximar la solucion de (6.32) en los nodos de la red (xj , tn ),
denotada por unj = u(xj , tn ) por valores numericos Ujn , j = 0, . . . , J, n = 0, . . . , N . Los ndices
j = 0, J representan nodos de la frontera, donde se imponen las condiciones de contorno, de
manera que
U0n = g0 (tn ), ULn = gL (tn ), n = 0, . . . , N.
Asmismo, el ndice n = 0 corresponde al instante inicial t0 = 0, donde la solucion esta dada por
u0 , de modo que
Uj0 = u0 (xj ), j = 0, . . . , J.
A partir de estos datos, la obtencion de la aproximacion numerica Ujn en los nodos internos viene
dada como solucion de una ecuacion en diferencias que aproxime a la ecuacion de (6.32) evaluada
en cada punto (xj , tn ) de la red. De esta manera, se calculan aproximaciones a las temperaturas
en los puntos xj seleccionados en la malla y en diferentes niveles de tiempo tn . La ecuaci
on en
diferencias depende de la aproximacion a las derivadas que elijamos.
Como ilustracion, consideraremos aqu diferentes casos de la familia de metodos conocida
como los -metodos [5]
!
n+1
n+1
n 2U n + U n
Ujn+1 Ujn
Uj+1
2Ujn+1 + Uj1
Uj+1
j
j1
= D
+ (1 )
k
h2
h2
+fjn+1 + (1 )fjn , j = 1, . . . , J 1, n = 0, . . . , N 1,
(6.35)
donde fjn = f (xj , tn ). La formula (6.35) se obtiene al aproximar la derivada temporal en (xj , tn )
por un cociente progresivo y la derivada segunda espacial por un promedio de diferencias centrales
en los puntos (xj , tn+1 ) y (xj , tn ), con [0, 1] un parametro. De la familia general, destacaremos
tres metodos.
M
etodo explcito. Corresponde al valor = 0 y la ecuacion en diferencias resultante tiene la
forma
n
n
Ujn+1 Ujn
Uj+1 2Ujn + Uj1
=D
+ fjn , j = 1, . . . , J 1, n = 0, . . . , N 1,
k
h2
o, alternativamente
n
n
Ujn+1 = rUj+1
+ (1 2r)Ujn + rUj1
+ kfjn , j = 1, . . . , J 1, n = 0, . . . , N 1,
(6.36)
n+1
v n
j1
j+1
M
etodo implcito
Corresponde al valor = 1 de la familia (6.35), es decir
!
n+1
n+1
Uj+1
2Ujn+1 + Uj1
Ujn+1 Ujn
=D
+ fjn , j = 1, . . . , J 1, n = 0, . . . , N 1,
k
h2
o, alternativamente
n+1
n+1
rUj1
+ (1 + 2r)Ujn+1 rUj+1
= Ujn + kfjn , j = 1, . . . , J 1, n = 0, . . . , N 1,
(6.37)
j1
v n+1
n
j+1
j
Figura 6.7: Metodo implcito
puede comprobarse que la matriz que resulta es estrictamente diagonalmente dominante, por
lo que queda asegurada la resolubilidad del sistema. Puede, ademas, utilizarse le algoritmo de
Thomas para su implementacion.
Pseudoc
odigo del m
etodo implcito.
Datos de entrada:
L: longitud del intervalo.
N : tama
no de la discretizacion temporal.
J: tama
no de la discretizacion espacial.
r: tama
no del cociente Dk/h2 .
U u0 . dato inicial.
f, g: terminos fuente y de frontera.
Inicializacion
h 1/J, k rh2 , T N k
Desde n = 1 hasta N
g0n g0 (tn ); gLn gL (tn )
Desde j = 1 hasta J 1
fjn f (xj , tn )
Vector del algoritmo de Thomas
p1 0
Desde j = 1 hasta J 1
pj+1 r/(rpj + (1 + 2r))
Resolucion del sistema
Desde n = 1 hasta N
U0 g0n ; UJ gLn
q1 U0 ;
Desde j = 1 hasta J 1
n + rq )/(rp + (1 + 2r))
qj+1 (U (j + 1) + kfj+1
j
j
Desde j = J 1 hasta 1
n +q
Uj pj Uj+1
j
M
etodo de Crank-Nicolson
El metodo correspondiente a = 1/2, o metodo de Crank-Nicolson, es de la forma
r n+1
r n+1
r n
r n
k n
Uj1
+ (1 + r)Ujn+1 Uj+1
=
Uj1 + (1 r)Ujn + Uj+1
+
fj + fjn+1 ,
2
2
2
2
2
j = 1, . . . , J 1, n = 0, . . . , N 1,
(6.38)
198
v n+1
v n
j1
j+1
El metodo tiene tambien un caracter implcito y puede implementarse a traves del algoritmo
de Thomas. Estas dos propiedades son, en realidad, retenidas por cualquier esquema de la familia
(6.35) con 6= 0.
Pseudoc
odigo del m
etodo Crank-Nicolson.
Datos de entrada:
L: longitud del intervalo.
N : tama
no de la discretizacion temporal.
J: tama
no de la discretizacion espacial.
r: tama
no del cociente Dk/h2 .
U u0 . dato inicial.
f, g: terminos fuente y de frontera.
Inicializacion
h 1/J, k rh2 , T N k
a r/2; b 1 + r; c 1 r
Desde n = 1 hasta N
g0n g0 (tn ); gLn gL (tn )
Desde j = 1 hasta J 1
fjn f (xj , tn )
Vector del algoritmo de Thomas
p1 0
Desde j = 1 hasta J 1
pj+1 a/(apj + b)
Resolucion del sistema
Desde n = 1 hasta N
U0 g0n ; UJ gLn
q1 U0 ;
Desde j = 1 hasta J 1
dj aUj1 + cUj aUj+1 + k2 (fjn + fjn+1
qj+1 (dj aqj )/(apj + b)
Desde j = J 1 hasta 1
n +q
Uj pj Uj+1
j
199
6.4.2.
An
alisis de convergencia
max ||U n un || = 0,
h,k0 0nkT
j = 1, . . . , J 1, n = 0, . . . , N 1,
con r = Dk/h2 . Cada componente de n se estima mediante desarrollos de Taylor en puntos
convenientes. Ello exige suficiente regularidad de la solucion teorica de (6.32). Se dice que el
metodo (6.35) es consistente si, para u suficientemente regular,
lm
max || n || = 0.
h,k0 0nkT
0nkT
k, h 0.
1
1
= ( )kuxxt h2 uxxxx
2
12
1 2
1 2
+ k uttt k uxxtt
24
8
1 1
2
+ ( )kh2 uxxxxt h4 uxxxxxx + ,
12 2
6!
+(1 2r(1
))enj
r(1
)enj+1 , j
(6.39)
= 1, . . . , J 1, n = 0, . . . , N 1,
0jJ
0jJ
entonces la hipotesis de estabilidad asegura que los coeficientes de (6.39) son no negativos y se
tiene
(1 + 2r)E n+1 2rE n+1 + E n + kT n+1/2 ,
de donde
E n+1 E n + kT n+1/2 .
Entonces, puesto que E 0 = 0,
n
E k
n1
X
m=0
u(0, t) = 0,
u(1, t) = 0,
(6.40)
t 0,
cuya solucion teorica es u(x, t) = et sin x. Las siguientes tablas corresponden a los errores en
norma dos de los tres esquemas: explcito (FTCS), implcito (BTCS) y Crank-Nicolson, aplicados
201
al problema de valores iniciales y de contorno (6.40). En todos los casos se han ejecutado los
programas hasta un tiempo final tf = 1 y las longitudes de paso se han dividido sucesivamente
por dos.
Comenzamos por el esquema explcito (FTCS). En la segunda columna de la tabla se han
calculado los errores en tiempo 1 y con = k/h2 2 constante e igual a 4/ 2 . Observemos que
entonces, si dividimos h por dos, el paso en tiempo k se divide por cuatro (o lo que es lo mismo:
si J se multiplica por dos, N ha de multiplicarse por cuatro). La tercera columna corresponde al
mismo esquema, pero con = 8/ 2 . Las diferencias son evidentes. Mientras en la primera columna
el metodo es convergente con orden dos (como es constante, k = O(h2 )), en la segunda, el valor
de rompe el esquema, haciendole perder su estabilidad.
||U u||2
M
etodo FTCS
0.1
0.05
0.025
0.0125
||U u||2
= 4/ 2
= 8/ 2
3.1060E-03
7.6843E-04
1.9161E-04
4.7872E-05
4.5265
8.6552E+52
1.5020E+263
En las tablas siguientes se tienen los errores para el metodo implcito (BTCS) y el metodo
Crank-Nicolson (CN). En ambos, = 4/ 2 . Los dos metodos son convergentes de orden dos y el
metodo Crank-Nicolson parece mas eficiente que el implcito, pues sus errores son menores.
M
etodo Crank-Nicolson
202
||U u||2
0.1
0.05
0.0025
0.0125
7.2170E-03
1.8276E-03
4.5838E-04
1.1469E-04
||U u||2
0.1
0.05
0.025
0.0125
2.1071E-03
5.3282E-04
1.3359E-04
3.3421E-05
= 4/ 2
= 4/ 2
6.4.3.
Problemas hiperb
olicos
Como modelo de una EDP hip`erbolica, consideramos un problema de valor inicial y contorno
asociado a la ecuacion de ondas
utt (x, t) = c2 uxx (x, t) + f (x, t),
0 < x < L, 0 t T,
0 < x < L,
u(0, t) = g0 (t),
0 t T.
u(L, t) = gL (t),
(6.41)
un+1
2unj + ujn1
j
+ O(k 2 ), k 0,
k2
unj+1 2unj + unj1
+ O(h2 ), h 0,
h2
n 2U n + U n
Uj+1
j
j1
=c
+ fjn , j = 1, . . . , J 1, n = 1, . . . , N 1, (6.42)
h2
2
k2
utt (x, 0) + O(k 3 )
2
k2 2
c uxx (x, 0) + f (x, 0) + O(k 3 )
2
k 2 2 00
k2
= u0 (x) + kv0 (x) + c u0 (x) + f (x, 0) + O(k 3 ),
2
2
siempre que u0 sea dos veces derivable. Este desarrollo proporciona la aproximacion
= u0 (x) + kv0 (x) +
k2 0
c2 k 2
0
0
0
U
2U
+
U
+ fj , j = 1, . . . , J 1,
j+1
j
j1
2h2
2
203
(6.43)
Uj1 = kVj0 +
Ujn+1
con = k/h.
Convergencia
Teniendo presente las formulas de aproximacion a las derivadas que hemos utilizado, puede
comprobarse que el esquema (6.44) es consistente de orden dos tanto en tiempo como en espacio.
Por otro lado, la ecuacion para los errores tiene la forma
en+1
= en1
+ c2 2 enj1 + 2(1 c2 2 )enj + c2 2 enj+1 + k 2 jn+1 ,
j
j
j = 1, . . . , J 1, n = 1, . . . , N 1,
con jn+1 el error de truncacion del esquema, obtenido de la manera usual, al sustituir en la
formula la solucion numerica por la teorica y dividir por k 2 . Por otra parte, la condici
on de
estabilidad del metodo es 0 < c 1. Si se satisface esta condicion y existen y son continuas las
derivadas utttt , uxxxx de la solucion teorica, el metodo es convergente de orden dos en tiempo y
en espacio [7, 9].
6.5.
Ejercicios
Ejercicio 1. Describe las formulas del metodo de Taylor de orden 2 para la ecuacion diferencial
dy/dt = f (t, y) = 1 + y 2 .
Ejercicio 2. La tabla 6.4 muestra, para diferentes longitudes de paso, el error en t = 3 de la
aproximacion proporcionada por un metodo de Taylor a la solucion del problema y 0 = (t y)/2
en [0, 3] con y(0) = 1.
h
1
0.5
0.25
0.125
error
7.9550E-04
5.0321E-05
3.3214E-06
2.0131E-07
h
(2k1 + 3k2 + 4k3 ),
9
204
donde
k1 f (t, y(t)), k2 f (t +
h
hk1
3h
3hk2
, y(t) +
), k3 f (t +
, y(t) +
).
2
2
4
4
Ejercicio 4. Comprueba que cuando se aplica el metodo de Euler para integrar la ecuaci
on
0
x (t) = x(t), la formula para avanzar un paso en la solucion es Yn+1 = (1 + h)Yn .
Ejercicio 5. Determina las condiciones que deben cumplir los coeficientes de un metodo RK
explcito de tres etapas para tener orden tres. Entre todos ellos, encuentra el que verifica c2 = c3
y b2 = b3 (conocido como metodo de Nystrom de tercer orden).
Ejercicio 6. Determina las condiciones que deben cumplir los coeficientes de un metodo RK
explcito de cuatro etapas para tener orden cuatro.
Ejercicio 7. Determina las condiciones para que un metodo explcito de dos etapas
k1 = f (tn + c1 h, yn )
k2 = f (tn + c2 h, yn + ha21 k1 )
yn+1 = yn + h(b1 k1 + b2 k2 ),
tenga orden dos, supuesto que no se verifican las condiciones c1 = 0 y c2 = a21 .
Ejercicio 8. Construye todos los metodos de orden dos de la forma
0
c2
c3
c2
0
0
c3
0
1/2
0
0
1/6
1/2
0
1/3
1
1/3
1/6
error
5.0110E-04
3.1040E-05
1.9401E-06
1.2120E-07
(iii) Comprueba que cuando se usa este metodo para resolver el problema y 0 (t) = f (t) en [a, b]
con y(a) = 0, el resultado es
N 1
hX
y(b)
f (tk + 4f (tk+1/2 ) + f (tk+1 ) ,
6
k=0
1/2
-1
1/6
2
2/3
.
1/6
y(0) = 1, y(/2) = 3.
Ejercicio 12. Muestra como se puede usar el metodo de tiro para resolver un problema de la
forma
y 00 (t) + p(t)y 0 (t) + q(t)y(t) = r(t),
0
ha y(a) ka y (a) = ,
a t b,
0
hb y(b) + kb y (b) = .
y(0) = 1, y(/6) = 5,
y(0) = 1, y 0 (0) = z,
a, b > 0,
(6.45)
y el esquema de aproximacion
Ujn+1 Ujn
k
+a
n Un
Uj+1
j1
2h
= b
n 2U n + U n
Uj+1
j
j1
.
h2
(6.46)
(i) Comprueba que el esquema (6.46) es consistente de orden uno en tiempo y orden dos en
espacio.
(ii) Comprueba que el esquema (6.46) puede reescribirse
n
n
Ujn+1 = (1 2b)Ujn + b(1 )Uj+1
+ b(1 + )Uj1
.
Ejercicio 16. M
etodos pseudoimplcitos. Considera el problema
ut = Duxx , a x b, 0 t T, c > 0,
a x b,
u(x, 0) = u0 (x),
u(a, t) = ga (t),
(6.47)
= D
= D
n U n + U n+1 U n+1
Uj+1
j
j1
j
h2
n+1
n Un
Uj+1
Ujn+1 + Uj1
j
h2
, 1 j J 1, n 0.
, 1 j J 1, n 0.
En el nivel de tiempo n = 0, se toma la condicion inicial exacta, y en los nodos frontera la soluci
on
numerica coincide con la teorica en todos los niveles de tiempo.
(i) Comprueba que, a pesar de ser aparentemente implcito, el esquema puede implementarse de
forma explcita.
(ii) Analiza la consistencia de los dos esquemas.
(iii) Analiza su estabilidad en norma del supremo cuando = ck/h2 es constante.
207
Referencias
208