Métodos Numéricos en Ingeniería
Métodos Numéricos en Ingeniería
2
Índice general
1
Métodos Numéricos Facultad de Ingenierı́a – IMERL
3. Interpolación polinomial 63
3.1. Introducción. ¿Por qué interpolar? . . . . . . . . . . . . . . . . . . . . . . . 63
3.2. Polinomio interpolante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.2.1. Forma de Vandermonde . . . . . . . . . . . . . . . . . . . . . . . . 66
3.2.2. Forma de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.2.3. Forma de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.3. Interpolación de funciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.3.1. Error de interpolación . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.3.2. Fenómeno de Runge . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.4. Interpolación a trozos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.4.1. Interpolación lineal a trozos . . . . . . . . . . . . . . . . . . . . . . 81
3.4.2. Interpolación cúbica a trozos . . . . . . . . . . . . . . . . . . . . . . 84
4. Ecuaciones no lineales 93
4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.2. Método de bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.3. Método de la regla falsa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.4. Método de la secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.5. Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.6. Resumen: una muy breve comparación . . . . . . . . . . . . . . . . . . . . 106
4.7. Métodos iterativos generales . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.7.1. Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.7.2. Orden y velocidad de convergencia . . . . . . . . . . . . . . . . . . 111
2
Métodos Numéricos Facultad de Ingenierı́a – IMERL
3
Métodos Numéricos Facultad de Ingenierı́a – IMERL
4
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Importante
Este es un material preparado para el curso de Métodos Numéricos. La primera
versión fue para el curso 2023, y para esta edición hemos hecho algunos ajustes: mo-
dificamos contendios de varias secciones y capı́tulos, agregamos un capı́tulo acerca
de integración numérica y otro sobre descomposición en valores singulares.
Si encontrás errores o tenés comentarios o sugerencias sobre estas notas,
te agradecemos escribir en el foro de Sugerencias y erratas en las notas
teóricas la página del curso.
Las notas están en su segundo semestre en uso. Varios estudiantes hicieron suge-
rencias y comentarios, encontraron errores de tipeo o inconsistencias, y con eso
ayudaron a mejorar las notas. De todos modos, es inevitable que sigan habiendo im-
perfectos; en la medida que hagamos correcciones sobre lo que ya está escrito, en el
texto van a estar en rojo y vamos a dejarlas detalladas en la lista de actualizaciones.
Notación
Pautamos aquı́ muy brevemente alguna notación que utilizamos en las notas del curso.
Cuando escribamos un sı́mbolo de igualdad con dos puntos por delante, esto denota
que se trata de una definición. Por ejemplo, si ponemos a := b + c (o b + c =: a)
estamos definiendo la variable a como la suma de b y c; en cambio, si ponemos
a = b + c (o b + c = a) esta relación se deduce de las definiciones de a, b, c.
Hacemos un uso un tanto liberal del sı́mbolo “≈”, que indica que estamos haciendo
una aproximación, ya sea de un valor numérico o de una función. Por ejemplo,
perfectamente podrı́amos escribir
t2
π ≈ 3, π ≈ 3,1416, o cos(t) ≈ 1 − para t ≈ 0.
2
Escribimos los números reales usando la fuente normal, y los vectores en negrita;
en general, usaremos minúsculas para todos ellos. Ası́, tendremos x ∈ R y x ∈ Rn .
En cambio, escribimos las matrices usando mayúsculas con fuente normal: A ∈
Mm×n (R). Al conjunto de matrices cuadradas n × n lo denotamos simplemente
como Mn (R).
5
Métodos Numéricos Facultad de Ingenierı́a – IMERL
f (x)
lı́m = 0.
x→a g(x)
6
Métodos Numéricos Facultad de Ingenierı́a – IMERL
7
Métodos Numéricos Facultad de Ingenierı́a – IMERL
8
Capı́tulo 1
9
Métodos Numéricos Facultad de Ingenierı́a – IMERL
estrategia de aproximación para estimar sus soluciones, lo que requiere el diseño de algo-
ritmos. Finalmente, seguramente realicemos un programa computacional que implemente
nuestro algoritmo y nos devuelva una solución computacional a nuestro problema.
Esquemáticamente, tenemos los siguientes pasos.
Para el modelado, seguramente utilicemos las leyes que nos da la mecánica new-
toniana. ¿Debemos incluir el efecto del rozamiento del aire? No hacerlo conduce a
un modelo muy simple (caı́da libre), pero que no serı́a apropiado si, por ejemplo,
queremos determinar la velocidad de caı́da de una hoja de un árbol. En algunos
casos, como en la caı́da de una gota de lluvia, puede ser relevante que el modelo in-
cluya que la masa del objeto pueda cambiar con el tiempo (porque la gota podrı́a ir
absorbiendo pequeñas gotitas del aire). Un buen modelo debe ser lo suficientemente
rico como para capturar los fenómenos que nos interesan con la precisión que nos
interesa, pero también lo suficientemente simple como para ser tratable.
Una vez determinadas las ecuaciones modelando la caı́da del objeto, éstas involucran
ciertos parámetros que debemos medir o estimar, como por ejemplo la altura
y velocidad inicial de la caı́da, o la constante gravitatoria. En caso de incluir el
rozamiento del aire, también podrı́an ser relevantes la forma y masa del objeto y
el coeficiente de rozamiento. Para una gota de lluvia cuya masa va variando, es
necesario incluir en el modelo la humedad del aire que la rodea, que también debe
ser determinada.
Nos encontramos ahora ante un sistema de ecuaciones (en este ejemplo, diferen-
ciales) que debemos resolver. En algún caso es posible que podamos resolverlas
analı́ticamente, pero en general necesitaremos una estrategia para aproximar las
soluciones. Esto es lo que se conoce como una discretización; hacia el final de es-
te curso aprenderemos métodos para tratar con ecuaciones diferenciales ordinarias.
Podemos aproximar las ecuaciones diferenciales involucradas mediante un sistema
de ecuaciones algebraicas si, por ejemplo, reemplazamos las derivadas por cocientes
incrementales,
f (x + h) − f (x)
f 0 (x) ≈ con h > 0 “chico”.
h
Finalmente, debemos resolver nuestro sistema de ecuaciones algebraicas. El tamaño
del sistema implica que seguramente sea inviable resolverlo “a mano”, y buscare-
mos implementar un programa computacional para realizar la tarea. La salida
10
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(output) de este programa deben ser las cantidades que deseemos estimar, como por
ejemplo la velocidad terminal de caı́da del objeto.
Este curso se concentra en los dos últimos puntos discutidos en el Ejemplo 1.1.1: el di-
seño, implementación y análisis de métodos para la solución de problemas de ingenierı́a.
Asumimos que tenemos dadas las ecuaciones o los problemas que queremos resolver. Bus-
camos desarrollar herramientas que nos permitan tratar computacionalmente con estos
modelos, comparar diferentes métodos entre sı́, entender en qué casos uno puede ser mejor
que otro, estimar sus costos computacionales, e introducir nociones básicas del área de la
matemática llamada Análisis Numérico.
1.1.1. Errores
Al resolver problemas reales mediante modelos computacionales, debemos lidiar con di-
versos errores. En cada una de las etapas de las descritas en el Ejemplo 1.1.1 se cometen
errores. En referencia a los cuatro puntos mencionados en ese ejemplo, podemos tener:
Además de estos errores, no podemos descartar los errores humanos (en cualquiera de las
etapas) y los errores de instrumentos (como fallas en el equipamiento para medir, defectos
de hardware). En este curso no vamos a tener en cuenta estos últimos factores.
En este capı́tulo nos centramos en los errores asociados a la representación de punto flo-
tante y de operaciones entre números representados en punto flotante. En el resto del
11
Métodos Numéricos Facultad de Ingenierı́a – IMERL
curso, nos vamos a enfocar en errores de aproximación para diversos métodos, pero en la
Sección 1.2 vamos a mostrar con algunos ejemplos a qué nos referimos con la expresión
“errores de aproximación”. La Sección 1.3 define los conceptos de error absoluto y relativo
y nos da algunas pautas sobre cómo los errores se pueden propagar al aplicar funciones u
operaciones aritméticas elementales. En la Sección 1.4 introducimos el sistema de nume-
ración que usan las computadoras, vemos algunas constantes especiales y consideramos
su capacidad de representar números reales. Finalmente, la Sección 1.5 se enfoca en cómo
los distintos tipos de error pueden aparecer y propagarse en un algoritmo.
está bien definida (el integrando es monótono creciente en el intervalo [0, 1]) pero el valor
de I no puede escribirse usando una cantidad finita de funciones elementales: no conocemos
una primitiva de esen x . Podemos considerar una partición uniforme del intervalo [0, 1] con
n + 1 puntos, Pn = {a0 , a1 , . . . , an }, esto es, ai = i/n con i = 0, . . . , n, y computar las
sumas inferiores y superiores asociadas a Pn ,
n−1 n−1
X
sen x
X 1 sen( ni )
S∗ (f, Pn ) = (ai+1 − ai ) ı́nf e = e ,
i=0
x∈[ai ,ai+1 ]
i=0
n
n−1 n−1
X X 1 sen( i+1
n ).
S ∗ (f, Pn ) = (ai+1 − ai ) sup esen x = e
i=0 x∈[ai ,ai+1 ] i=0
n
12
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por ejemplo, usando una partición P10 con n = 10 intervalos, esto nos conduce a
Esto nos da cotas bastante groseras para el valor de I: haciendo el promedio entre nuestra
cota superior y nuestra cota inferior, deducimos que I ≈ 1,632 ± 0,066. Si esta precisión
fuese suficiente para los cálculos que necesitamos, podrı́amos quedarnos satisfechos. Si no
lo fuese, podrı́amos tomar mayor cantidad de puntos.
Ejercicio 1.1. Implementar un código de Octave que permita reproducir los resultados
de la tabla de arriba.
sen(k) (0) −1 k
sen(10−1 ) = sen(0) + cos(0)10−1 + . . . + (10 ) + rk (10−1 ),
k!
donde
sen(k+1) (θx ) −1 k+1
rk (10−1 ) = (10 ) , y θx ∈ [0, 10−1 ].
(k + 1)!
13
Métodos Numéricos Facultad de Ingenierı́a – IMERL
14
Métodos Numéricos Facultad de Ingenierı́a – IMERL
n an a2n − 7
0 3 2
8
1 3
= 2,6666̄ 0.1̄
2 127
48
= 2,6458 3̄ ≈ 4,3 × 10−4
3 32257
12192
≈ 2,645751312335958 ≈ 6,7 × 10−9
1.3. Errores
La noción de error se refiere a una discrepancia entre lo que computamos y lo que queremos
computar. En esta sección introducimos la diferencia entre error absoluto y relativo y
estimamos cómo una discrepancia en las entradas afecta a la salida al evaluar una función
o al realizar operaciones aritméticas elementales (suma, resta, producto, división).
ex := x̄ − x.
15
Métodos Numéricos Facultad de Ingenierı́a – IMERL
4
Definición 1.3.2 (error relativo). Sean x ∈ R un número real dado con x 6= 0 y x̄ ∈ R
una aproximación de x. Definimos el error relativo como el error absoluto normalizado
por x,
x̄ − x ex
εx := = .
x x
4
Observación 1.3.1 (errores con signo). Las definiciones anteriores implican que los errores
absoluto y relativo tienen signo: pueden ser positivos o negativos. Elegimos esta convención
porque nos permite relacionar las cantidades x y x̄ mediante las sencillas identidades
x̄ = x + ex y
x̄ = (1 + εx ) x. (1.2)
En general, en las aplicaciones nos interesa la magnitud de los errores y no su signo. En
esos casos, tenemos que estimar o acotar |ex | o |εx |. 4
16
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(x̄−x)
Al ser x
= εx , deducimos que
f 0 (x)x
εy ≈ εx ,
f (x)
lo que significa que el error relativo en y = f (x) se relaciona con el error relativo en x
mediante un cierto factor. Esto motiva la siguiente definición.
Definición 1.3.3 (número de condición). Llamamos número de condición de f en
el punto x al número
f 0 (x)x
κf (x) := . (1.3)
f (x)
Vagamente,
si κf (x) ∼ 1, diremos que el problema de evaluar f en x está bien condicionado;
κf (x) 1, diremos que el problema de evaluar f en x está mal condicionado.
4
Remarcamos que el hecho de que el problema de evaluar f en x esté bien o mal con-
dicionado solamente tiene que ver con la función f . No estamos asumiendo nada sobre
cómo estamos computando f . En matemática, los números de condición dan una pauta
de la sensibilidad en la salida de una operación (como evaluar una función) respecto a sus
entradas. Al trabajar con sistemas de ecuaciones lineales, vamos a retomar este concepto
y dar una definición análoga a (1.3).
17
Métodos Numéricos Facultad de Ingenierı́a – IMERL
División. Suponiendo que ni y ni ȳ son nulos, dividimos las dos igualdades en (1.4) y
1
usamos el equivalente1 1+t ≈ 1 − t para obtener
x̄ x (1 + εx ) x x
= ≈ (1 + εx )(1 − εy ) ≈ (1 + εx − εy ).
ȳ y (1 + εy ) y y
εx/y ≈ εx − εy .
En consecuencia, en el peor de los casos, el error relativo en la división entre dos números
es comparable en magnitud al peor de los errores relativos en ellos. Recordamos que los
errores relativos son tomados con signo, por lo que el signo negativo en εy no tiene ningún
significado especial.
x̄ + ȳ − (x + y) x y
εx+y = = εx + εy .
x+y x+y x+y
x y
Aquı́, vemos que aparecen los factores x+y
e x+y
multiplicando a εx y εy , respectivamente.
Distinguimos tres casos.
Si xy > 0 (es decir, si x e y tienen igual signo), entonces los factores que multiplican
a εx y εy son ambos positivos, suman 1, y podemos acotar
x y
|εx+y | ≤ máx{|εx |, |εy |} + máx{|εx |, |εy |} = máx{|εx |, |εy |}.
x+y x+y
Si xy < 0 y |x| |y| (o |x| |y|), entonces los factores que multiplican a εx y εy
no son mucho mayores que 1 en valor absoluto y tenemos
x y
|εx+y | ≤ |εx | + |εy | . |εx | + |εy |.
x+y x+y
En este caso tenemos que, en el peor de los casos, el error relativo en la suma x + y
es comparable en magnitud con el peor de los errores relativos en x e y.
1 1
Este equivalente no es otra cosa que el polinomio de Taylor de orden 1 de la función f (t) = 1+t
alrededor del origen.
18
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Vamos a ilustrar por qué se usa el tan dramático término “cancelación catastrófica” con
un ejemplo.
Ejemplo 1.3.1 (cancelación catastrófica). Buscamos las raı́ces del polinomio x2 − 56x +
1 = 0. Es fácil mostrar que estas raı́ces son
√ √
r1 = 28 + 783 ≈ 55,982; r2 = 28 − 783 ≈ 0,017863.
Supongamos que √ solamente tenemos 5 cifras de precisión. Con ellas, 28 se computa
exactamente y 783 se aproxima como 27,982. Tenemos errores relativos ε28 = 0 y
ε√783 ≈ −4,90 × 10−6 .
√
La raı́z r1 se puede calcular sin error considerable (e incluso mejor que el de 783) pues
estamos sumando números positivos,
r1 = 28 + 27,982 = 55,982 εr1 ≈ −2,45 × 10−6 .
Sin embargo, el cálculo de la raı́z r2 implica una resta de dos valores muy próximos,
r2 = 28 − 27,982 = 0,018, εr2 ≈ 7,7 × 10−3 .
Observamos que el error relativo en r2 aumentó en magnitud en un factor de aproxima-
damente 1500 respecto al error relativo en las entradas que utilizamos. Este aumento de
magnitud
√ puede causar que un error tolerable en las entradas (nuestras aproximaciones
de 28 y 782) dé lugar a un error no tolerable en la aproximación de r2 : esto es lo que se
llama una cancelación catastrófica.
Para evitar caer en una cancelación catastrófica, es importante evitar computar r2 en la
forma en que lo hicimos. Una forma de lograrlo es reescribir el polinomio original como
x2 − 56x + 1 = (x − r1 )(x − r2 ) = x2 − (r1 + r2 )x + r1 r2 .
Podemos determinar r2 igualando los coeficientes término a término y usando el valor ya
calculado para r1 . Si utilizamos el término independiente se tiene
1 1
r1 r2 = 1 ⇔ r2 = = ≈ 0,017863, εr2 ≈ 8,9 × 10−6 .
r1 55,982
Esta es una aproximación que tiene un error√relativo del mismo orden de magnitud que
el que tenemos en nuestra aproximación de 783. En cambio, observar que si en lugar
del término independiente quisiéramos determinar r2 usando el término en x, entonces
llegarı́amos al mismo mal resultado que antes. 4
Del ejemplo anterior, sacamos una importante moraleja: cuando nos enfrentamos a
problemas que implican restar números de igual magnitud, es importante bus-
car reescribir las fórmulas para evitar cancelaciones catastróficas.
19
Métodos Numéricos Facultad de Ingenierı́a – IMERL
¿Por qué ocurre esto? La máquina realiza las operaciones en el orden en que se le presentan,
por lo que primero suma 1 + 10−16 y al resultado que obtiene le resta 1. El hecho de que el
resultado final sea igual a 0 quiere decir que, para la máquina, los números 1 + 10−16
y 1 son indistinguibles.
En general, al almacenar un número x podremos guardar algún número próximo a éste
de entre el conjunto de números que puede representar exactamente la máquina. A conti-
nuación estudiamos este conjunto y algunas de sus propiedades básicas, luego analizamos
qué efectos puede tener este error de representación (o almacenamiento), y finalmente
marcamos algunos cuidados que debemos tener en consecuencia.
20
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ası́, los números de punto flotante son aquellos que se pueden escribir de la forma
x = ±(1 + f ) × 2e . (1.5)
Aquı́, la mantisa es un número 0 ≤ f < 1 tal que 252 f es un número entero, esto es,
j
f= , j ∈ {0, 1, . . . , 252 − 1}.
252
Notemos que hay exactamente 252 posibles valores para f , lo que se corresponde con el
hecho de que f sea representado por 52 bits.
En cambio, el exponente e es un número entero tal que −1022 ≤ e ≤ 1023. Observemos
que esta elección nos deja 2046 = 211 − 2 posibles valores para e; los dos valores restantes
para el exponente (e = −1023 y e = 1024) se reservan para ciertos números especiales.
Distribución de números
Analicemos con más detalle cómo se distribuyen los números de punto flotante. Para fijar
ideas, consideremos los posibles números de punto flotante que se corresponden con signo
positivo y e = 0 en (1.5). Al variar la mantisa f , tendremos un conjunto de 252 números,
252 − 1
1 2 3
1, 1 + 52 , 1 + 52 , 1 + 52 , . . . 1 + .
2 2 2 252
Este conjunto es el conjunto de números en el intervalo [1, 2) que pueden ser representados
exactamente por una máquina que use el estándar de precisión doble. Notamos que este
conjunto es un conjunto de 252 números equiespaciados, con una separación igual a 2−52 .
Ahora, fijemos un valor de e mayor en (1.5), por ejemplo e = 4. Al variar la mantisa
f , nuevamente tendremos un conjunto de 252 números equiespaciados pero ahora en el
intervalo [24 , 25 ) = [16, 32),
252 − 1
1 2 3
16, 1 + 52 × 16, 1 + × 16, 1 + 52 × 16, ... 1 + × 16 .
2 252 2 252
La separación entre elementos consecutivos de este conjunto es 16
252
= 2−48 , esto es, 16
veces mayor que la separación en el intervalo [1, 2).
Observamos, pues, que el conjunto de números de punto flotante es finito y la
distancia entre números consecutivos no es globalmente uniforme.
Redondeo y truncamiento.
Como hemos visto, los sistemas de representación pueden representar una cantidad finita
de números, por lo que al desear representar un número real x que quizás no esté incluido
en el mismo, la computadora deberá aproximarlo con un número fl(x) que sı́ lo esté. Para
esta aproximación existen dos métodos habitualmente usados:
21
Métodos Numéricos Facultad de Ingenierı́a – IMERL
lo que arroja el número x = 1 + 2−52 en (1.5). Deducimos, por lo tanto, que en precisión
doble
εM = 2−52 ≈ 2,2 × 10−16 .
En Octave y Matlab, el comando eps nos devuelve el valor de εM para precisión doble,
que es el formato que nuestras máquinas usan por defecto:
Como los sistemas de punto flotante incluyen una cantidad finita de números, necesaria-
mente debe haber uno máximo y uno mı́nimo (en valor absoluto) representables.
22
Métodos Numéricos Facultad de Ingenierı́a – IMERL
lo que da lugar a
Realmin = 1 × 2−1022 ≈ 2,2 × 10−308 .
En Octave y Matlab, los comandos realmax y realmin nos devuelven los valores de
Realmax y Realmin para precisión doble, respectivamente:
4
Recordamos que los exponentes e = −1023 y e = 1024 se reservan para ciertos números
especiales. Veamos algunos de estos casos.
Definición 1.4.3 (overflow e Inf). Si un cálculo en la máquina arroja un número mayor
que Realmax , decimos que se produce un overflow. El resultado de este cálculo es un
número de punto flotante llamado Inf, que se almacena con el exponente e = 1024 y la
mantisa f = 0. 4
Definición 1.4.4 (NaN). Si hacemos un cálculo en la máquina que no está definido en
el sistema real, como por ejemplo 0/0 o Inf − Inf, el resultado de nuestro cálculo es un
número de punto flotante llamado NaN2 , que se almacena con el exponente e = 1024 y
una mantisa f 6= 0. 4
23
Métodos Numéricos Facultad de Ingenierı́a – IMERL
obtenemos que ambos resultados son cero. Ahora, si usamos estas variables para calcular:
tenemos dos infinitos distintos. En particular, escribiendo el cero positivo como +0,0 y el
negativo como −0,0, en el estándar IEEE 754 se cumple
1 1
= +Inf, = −Inf.
+0,0 −0,0
Por defecto, al computar o escribir 0 en Octave o Matlab, obtenemos +0,0.
24
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Demostración. Sea x tal que 2−1022 ≤ |x| < 21024 ; sin perder generalidad, asumimos que
x es positivo. Consideremos los intervalos de la forma [2e , 2e+1 ) con e entero y tal que
−1022 ≤ e ≤ 1023. Como estos intervalos son disjuntos dos a dos y la unión de todos
ellos es [2−1022 , 21024 ), necesariamente x tiene que estar en uno de ellos. Esto es, existe un
único e entero con −1022 ≤ e ≤ 1023 y αx ∈ [0, 1) tales que
x = (1 + αx ) × 2e .
−52
Si αx ≤ 1 − 2 2 , entonces la representación de punto flotante de x usando redondeo es
un número fl(x) de la forma
fl(x) = (1 + fx ) × 2e .
Notemos que el exponente es el mismo que en x, y la mantisa es fx = k2−52 para algún
k = 0, . . . , 252 − 1. Si restamos, dividimos por x y simplificamos las exponenciales, nos
encontramos con que
fl(x) − x fx − αx
= .
x 1 + αx
Finalmente, tenemos que 1 + αx ≥ 1, y como estamos usando redondeo se cumple |αx −
−52
fx | ≤ 2 2 ; notar que como el espaciado entre mantisas es de 2−52 , la mantisa más cercana
a αx no puede estar a una distancia mayor que la mitad de este espaciado. Tomando valor
absoluto en la igualdad anterior y acotando, se llega a (1.6).
2−52
Por otra parte, si 1 − 2
< αx < 1, entonces la representación de punto flotante de x
usando redondeo es
fl(x) = 1 × 2e+1 .
Dejamos como ejercicio verificar que en ese caso también se cumple (1.6).
Observación 1.4.3 (cota inferior para |x|). La condición sobre |x| en la Proposición 1.4.1
es necesaria para evitar que fl(x) sea Inf o 0, aunque la cota inferior puede mejorarse si
se usan números subnormales. También observemos que esta proposición refiere al error
relativo al aproximar con redondeo un número en punto flotante; si nos interesa estimar
el error absoluto en esta aproximación, usando la Definición 1.3.1 tenemos
|x| εM
|ex | = |x εx | ≤
.
2
Por lo tanto, mientras nuestra mejor cota para el error relativo al aproximar un
número real es uniforme, nuestra mejor cota para el error absoluto depende
de la magnitud de dicho número real. 4
25
Métodos Numéricos Facultad de Ingenierı́a – IMERL
1.5.1. Advertencias
La aritmética en punto flotante no es ni asociativa ni distributiva: las operaciones se
hacen por etapas y en cada operación se aplica el redondeo correspondiente. Por lo tanto,
colocar paréntesis puede alterar el resultado en punto flotante incluso cuando no lo hace
en aritmética real. Recordando que en el formato de precisión doble εM = 2−52 ≈ 2,204 ×
10−16 , veamos un ejemplo.
¿Por qué obtenemos tres resultados diferentes3 ? En el primer caso, la operación entre
paréntesis da como resultado 1 (porque 1 + 10−16 se redondea a 1), al sumarle a ese 1
el número 10−16 se vuelve a obtener 1, y al restar 1 se obtiene 0. En el segundo caso,
primero se resuelve la suma entre paréntesis, cuyo resultado es 2.0000e-16. Al sumar 1
más dicho número el resultado es 1+eps (porque estamos usando redondeo), y al restar 1 se
obtiene como resultado final eps. Finalmente, en el tercer caso se resuelven las sumas entre
paréntesis en primer lugar: el resultado de la primera es 1 y el de la segunda es -1+0.5*eps
(ver el ejercicio a continuación). Al sumar ambas, el resultado final es 0.5*eps.
Ejercicio 1.4. Explicar por qué al computar 1e-16 - 1 en Octave se obtiene el resultado
-1+0.5*eps.
26
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Al término ε2M + κf (x) ε2M lo llamamos error inevitable al computar f (x). Es la mejor
cota que podemos tener para el error relativo en caso de que pudiésemos computar f
exactamente. Para valorar la calidad de un algoritmo, más que utilizar como referencia el
error en su salida, debemos comparar este error con el error inevitable. Si un algoritmo
para computar y = f (x) produce una salida ȳ con error relativo comparable con el error
inevitable, diremos que el algoritmo es “bueno”. Esta no es una definición matemática,
pero nos da una pauta de la calidad de un método dado. Veamos dos ejemplos.
Ejemplo 1.5.1 (algoritmo malo, algoritmo bueno). Queremos evaluar y = 1 − cos(x) en
x = 10−10 . Estimamos el número de condición de f (x) = 1 − cos(x) en el punto 10−10 (ver
(1.3)) usando que sen t ≈ t y 1 − cos t ≈ t2 /2 para t ≈ 0,
− sen(10−10 ) 10−10
κf (10−10 ) = ≈ 2.
1 − cos(10−10 )
Por (1.7), esto implica que el error inevitable al computar y = 1 − cos(10−10 ) es aproxima-
damente 32 εM ; en otras palabras, teóricamente deberı́a ser posible computar y sin cometer
un gran error relativo.
Para computar y, la estrategia que parece más obvia es la que harı́amos si tuviésemos
lápiz y papel: en primer lugar calcular y1 = cos(10−10 ), y luego computar y = 1 − y1 . Esto
se corresponde con componer f2 ◦ f1 (10−10 ), donde
Analicemos la mejor cota que podemos obtener para la magnitud del error relativo en
este algoritmo. Razonando como hicimos para llegar a (1.7), tenemos
εM
|εy | ≤ + κf2 (y1 ) |εy1 |,
2
y
εM εM
|εy1 | ≤ + κf1 (x) .
2 2
Aquı́ utilizamos que la magnitud del error relativo al aproximar x se puede acotar por ε2M ,
ya que x es la entrada de nuestro algoritmo. Combinando las dos expresiones, obtenemos
εM ε
M εM
|εy | ≤ + κf2 (y1 ) + κf1 (x) .
2 2 2
27
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por lo tanto, nos basta con estimar los números de condición κf1 (x) y κf2 (y1 ). Usando
nuevamente la definición (1.3) y las aproximaciones sen(10−10 ) ≈ 10−10 y 1 − y1 = 1 −
cos(10−10 ) ≈ 10−20 /2, tenemos
− sen(10−10 ) 10−10
κf1 (x) = ≈ 10−20 ,
cos(10−10 )
y1
κf2 (y1 ) = ≈ 2 × 1020 .
1 − y1
Sustituyendo, obtenemos la cota
εM ε
M εM
|εy | ≤ + 2 × 1020 + 10−20 ≈ εM × 1020 .
2 2 2
¡Esta cota es del orden de 1020 veces peor que el error inevitable! Esto nos da la pauta de
que este es un mal algoritmo: si lo implementamos computacionalmente, estamos intro-
duciendo cancelaciones catastróficas. Es muy sencillo verificarlo computacionalmente:
28
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ejercicio 1.5. Completar los detalles de la estimación del error para el segundo algoritmo
en el ejemplo anterior.
Ejemplo 1.5.2 (otro ejemplo de cancelaciones). Consideremos el polinomio
Es fácil verificar que p tiene una raı́z en x = 1. Corremos en Octave unas lı́neas de código
que nos permitan graficarlo para x ∈ (0,99, 1,01):
El resultado se muestra a la derecha en la Figura 1.1, y notamos que ahora sı́ se parece a
lo que esperamos del gráfico de un polinomio.
En teorı́a, las dos opciones que presentamos arriba nos permitirı́an graficar al polinomio
p en el intervalo (0,99, 1,01). Sin embargo, al evaluar p(x) con x cercano a 1, la primera
introduce cancelaciones que se pueden evitar si se factoriza a p. La segunda opción es
computacionalmente mucho más robusta que la primera. 4
Ejercicio 1.6. Justificar la última afirmación del ejemplo anterior: realizar un análisis de
error hacia adelante para computar p(x) con x cercano a 1 usando ambos métodos.
Observación 1.5.1 (análisis de error hacia atrás). El análisis de error hacia adelante mide
la discrepancia entre los valores computados y los valores exactos que queremos aproximar.
Un enfoque alternativo para estimar el error consiste en, dado un algoritmo, preguntarse
cuánta deberı́a ser la discrepancia entre la entrada que usamos y la entrada exacta para
que se produzca la salida aproximada que obtuvimos. Esta alternativa se suele llamar
análisis de error hacia atrás, pero no profundizaremos en ella. Para quien esté interesado,
referimos a [Hea02, Sección 1.2.5]. 4
29
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 1.1: Gráfico de p(x) usando (1.8) (izquierda) y usando (1.9) (derecha).
30
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Podrı́amos haber puesto fl(h) en vez de h en el denominador, pero como no vamos a tener
cancelaciones, la simplificación que hicimos no es relevante.
Estimemos la magnitud del error absoluto en la aproximación (1.11). Observamos que hay
dos fuentes de error:
Una forma de cuantificar esta observación es utilizar una desigualdad triangular para
separar el error absoluto como la suma de las componentes de truncamiento y de redondeo,
Error de truncamiento
Comenzamos por acotar el primer término del lado derecho en (1.12). Este término tiene
que ver exclusivamente con reemplazar la derivada por un cociente incremental. Para
analizarlo, basta con utilizar un desarrollo de Taylor de primer orden en el punto x,
usando la expresión de Lagrange para el resto,
h2
f (x + h) = f (x) + f 0 (x) h + f 00 (c) para algún c ∈ [x, x + h].
2
Con esta expresión, podemos despejar y estimar
h
|errortruncamiento | = |f 00 (c)| .
2
Como estamos asumiendo que f es de clase C 2 y c ∈ [x, x + h], entonces esperamos que
f 00 (c) ≈ f 00 (x). Obtenemos ası́ una estimación para el error de truncamiento,
h
|errortruncamiento | ≈ |f 00 (x)| . (1.13)
2
Deducimos que el error de truncamiento es de orden h: esperamos que decrezca linealmente
al decrecer el paso h.
31
Métodos Numéricos Facultad de Ingenierı́a – IMERL
con
εM εM
|εf (x+h) | ≤ , |εf (x) | ≤ .
2 2
Por lo tanto, podemos reemplazar fl(f (x+h)) y fl(f (x)) y utilizar la desigualdad triangular
y las cotas para los errores relativos para obtener
|f (x + h) εf (x+h) − f (x) εf (x) |
|errorredondeo | =
h
|f (x + h)| |εf (x+h) | + |f (x)| |εf (x) |
≤
h
|f (x + h)| + |f (x)| εM
≤ .
h 2
Como f es continua y f (x) 6= 0, podemos suponer f (x + h) ≈ f (x), y concluimos que
|f (x)| εM
|errorredondeo | . . (1.14)
h
A diferencia del error de truncamiento, nuestra cota para este error es inversamente pro-
porcional al paso h. Esto significa que un paso demasiado pequeño puede incrementar el
error relacionado a la representación en punto flotante.
Paso óptimo
Tenemos entonces que, por un lado, se requiere un paso pequeño para minimizar el error de
truncamiento, pero por otro lado un paso demasiado pequeño puede dar lugar a problemas
de punto flotante. Una forma de encontrar el paso óptimo es combinar (1.13) y (1.14),
sustituir en (1.12) para obtener
h |f (x)| εM
|ef 0 (x) | ≤ |f 00 (x)| + ,
2 h
y buscar el valor hopt > 0 en el que el lado derecho se hace lo menor posible. Esto se logra
de forma sencilla derivando el lado derecho de la desigualdad respecto a h y despejando.
Obtenemos s
2 |f (x)| εM
hopt = , (1.15)
|f 00 (x)|
lo que da lugar a una cota óptima para la magnitud del error absoluto,
p
|ef 0 (x),opt | = 2 |f (x)| |f 00 (x)| εM (1.16)
32
Métodos Numéricos Facultad de Ingenierı́a – IMERL
En la Figura 1.2 ilustramos cómo se comportan las dos componentes del error y su suma al
variar h. Esta sencilla figura es ilustrativa del comportamiento cualitativo de los errores de
truncamiento y de redondeo, pero en la práctica podemos ver algo un poco más interesante.
00
etrunc = |f (x)|
2
h
ered = |f (x)|
h
εM
Figura 1.2: Cota superior del error total. Permite hallar h óptimo.
Ejemplo 1.5.3 (estimación del paso óptimo). Apliquemos la discusión anterior a apro-
ximar la derivada de la función f (x) = ex en el punto x = 1. Como f 00 (1) = f (1) = e, las
expresiones (1.15) y (1.16) nos indican que
√
hopt = 2εM = 2−52/2 ≈ 2,1 × 10−8 , |ef 0 (x),opt | = hopt e ≈ 5,7 × 10−8 .
Implementamos en Octave un código que tome varios h en el intervalo [10−16 , 10−1 ], calcule
las aproximaciones (1.11) y estime el error sabiendo que f 0 (x) = e.
n = 1000;
derivadas = zeros(1,n);
h = 10.^(-linspace(1,16,n)); % toma n puntos entre 10^{-16} y 10^{-1}
f = @(x) exp(x);
x = 1;
for k = 1:n
derivadas(k) = (f(x+h(k))-f(x))/(h(k));
end
33
Métodos Numéricos Facultad de Ingenierı́a – IMERL
La salida de este código se muestra en la Figura 1.3. Observamos que, para h grande,
el error computado decrece linealmente al reducir h: esta es la región en la que domina
el error de truncamiento. En cambio, para h pequeño el error se comporta de manera
errática aunque muestra una tendencia a crecer al reducir h: en esta región domina el
error de redondeo. La transición entre una región y la otra se da cuando h (y el error) son
del orden de 10−8 , consistentemente con nuestras estimaciones de hopt y |ef 0 (x),opt |.
Ejercicio 1.7. En la Figura 1.3 se observa que para h ≈ 10−16 ≈ εM el error absoluto es
del orden de una unidad. Estimar el error de redondeo para h = 10−16 .
34
Capı́tulo 2
2.1. Introducción
Un sistema de m ecuaciones lineales y n incógnitas consiste en un conjunto de relaciones
algebraicas de la forma
Xn
aij xj = bi ∀ i = 1 . . . m,
j=1
Cuando m = n, decimos que el sistema es cuadrado; en ese caso, el sistema tiene solución
única si y sólo si det A 6= 0. En este capı́tulo trabajamos con sistemas de este tipo, a los
que llamamos compatibles determinados.
Asumiendo que el sistema Ax = b tiene solución única, nos planteamos cómo hallarla.
En general, los métodos de resolución de sistemas lineales pueden ser clasificados en las
siguientes dos categorı́as.
35
Métodos Numéricos Facultad de Ingenierı́a – IMERL
det(Ai )
xi = ,
det(A)
5 1 3 5
det(A) = −4, det(A1 ) = = −4, det(A2 ) = = −8,
−1 −1 1 −1
Definición 2.1.1 (f lop1 ). Denotaremos por f lop a una simple operación de punto flo-
tante, como suma, resta, producto o división. 4
En general, no nos interesa la cantidad exacta de f lops que requiere un método, sino el
orden de magnitud de esta cantidad. Retomando el Ejemplo 2.1.1, supongamos que usamos
la regla de Cramer para resolver un sistema con A ∈ Mn (R). Calcular el determinante
de una matriz n × n desarrollando por filas o columnas requiere O(n!) operaciones, y la
regla de Cramer requiere que calculemos n + 1 determinantes de matrices n × n y los
dividamos n veces. Deducimos que una aplicación ingenua de la regla de Cramer para
resolver un sistema de n × n tiene un costo computacional de O((n + 1)!) f lops. Esto es
pésimo computacionalmente. Incluso aunque sea posible calcular determinantes de forma
eficiente (ver Observación 2.2.2 más abajo), es poco común usar la regla de Cramer para
resolver sistemas grandes de ecuaciones lineales. No vamos a profundizar en este algoritmo.
En la Sección 2.2 vamos a analizar con detalle otro método directo ya conocido de cursos
anteriores: el algoritmo de escalerización gaussiana. Luego, la Sección 2.3 profundiza sobre
la descomposición LU , que está cercanamente relacionada con dicho algoritmo. Si bien la
1
La palabra f lop proviene del inglés floating point operations.
36
Métodos Numéricos Facultad de Ingenierı́a – IMERL
37
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Como estamos asumiendo que este sistema tiene solución única, necesariamente tiene que
ser det(A) 6= 0. Al ser A triangular inferior, es fácil verificar que
de modo que aii 6= 0 para todo i = 1, . . . , n. Para resolver este sistema, parece natural
comenzar por la primera ecuación en (2.1):
b1
a11 x1 = b1 ⇒ x1 = .
a11
Luego, podemos usar el valor de x1 que hallamos para continuar con la segunda ecuación
en (2.1),
b2 − a21 x1
a21 x1 + a22 x2 = b2 ⇒ x2 = .
a22
Podemos continuar sucesivamente: si conocemos x1 , . . . , xi−1 , entonces podemos despejar
xi de la i-ésima ecuación en (2.1) mediante
Pi−1
bi − j=1 aij xj
xi = .
aii
Este método para resolver sistemas triangulares inferiores se denomina sustitución hacia
adelante. En el Algoritmo 2.1, lo presentamos en formato de pseudo-código.
Algoritmo 2.1: Pseudo-código: sustitución hacia adelante
Datos: A = (aij ) ∈ Mn (R) triangular inferior, b ∈ Rn
Resultado: x ∈ Rn solución de Ax = b
x1 ← ab111 ;
para i = 2 : n hacer
xi ← a1ii (bi − i−1
P
j=1 aij xj );
fin
bn
xn = ,
ann
y para i = n − 1, . . . 1 calcular
Pn
bi − j=i+1 aij xj
xi = .
aii
38
Métodos Numéricos Facultad de Ingenierı́a – IMERL
39
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(k)
aik
ese fin, para cada i = k + 1, . . . , n definimos el coeficiente lik := (k) , al que llamamos
akk
multiplicador, y actualizamos la fila i-ésima de la matriz A(k) y la entrada i-ésima del
término independiente b(k) mediante
Pn−1
loop en k: k=1 [2(n − k)2 + 3(n − k)] f lops.
n−1 n−1
X (n − 1)n X (n − 1)n(2n − 1)
k= , k2 = ,
k=1
2 k=1
6
40
Métodos Numéricos Facultad de Ingenierı́a – IMERL
obtenemos:
n−1
X
2n2 − 4nk + 2k 2 + 3n − 3k
f lops(EG) =
k=1
n−1
X n−1
X
2
k2
= 2n + 3n (n − 1) − (4n + 3) k+2
k=1 k=1
2n3
2
(n − 1)n (n − 1)n(2n − 1)
= 2n + 3n (n − 1) − (4n + 3) +2 =O .
2 6 3
2n3 2n3
2
O + O(n ) = O f lops.
3 3
Desde el punto de vista de la aritmética real, esto es, si tuviésemos precisión infinita,
(k)
pivotear sólo es necesario cuando se llega a que akk = 0 para algún paso k. Cuando
(k)
trabajamos con aritmética de punto flotante, pivotear se vuelve importante cuando akk ≈
0. Ilustramos esta idea con un ejemplo.
41
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Vale la pena comparar la solución exacta del problema, [0, −1, 1]t , con la que obtuvimos,
[−0,28, −1,4, 0,9993]t : comenzamos con un muy pequeño error en x3 , que se propagó
peligrosamente al realizar la sustitución hacia atrás. 4
42
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Una vez hallado el nuevo pivote, se intercambia la fila en la que está con la fila k-ésima
en la matriz ampliada y se continúa con la escalerización gaussiana.
Algoritmo 2.4: Pseudo-código: EG con pivoteo parcial
Datos: A = (aij ) ∈ Mn (R), b ∈ Rn
Resultado: A = (aij ) ∈ Mn (R) triangular superior, L = (lik ) ∈ Mn (R) triangular
inferior y que tiene los multiplicadores del algoritmo de escalerización
gaussiana sin pivoteo
para k = 1 : n − 1 hacer
p ← arg máxr∈{k,...,n} |ark |
intercambiar(bk , bp );
para j = 1 : n hacer
intercambiar(akj , apj );
fin
para i = k + 1 : n hacer
l(i, k) ← a(i, k)/a(k, k);
a(i, k) ← 0;
para j = k + 1 : n hacer
a(i, j) ← a(i, j) − l(i, k) ∗ a(k, j);
fin
bi ← bi − l(i, k) ∗ bk ;
fin
fin
43
Métodos Numéricos Facultad de Ingenierı́a – IMERL
2.3. Descomposición LU
En la Sección 2.2 analizamos el algoritmo de escalerización gaussiana (sin pivoteo y
con pivoteo parcial) para la resolución de sistemas de ecuaciones lineales Ax = b con
A ∈ Mn (R). Una consecuencia directa de este algoritmo es que nos da una forma intere-
sante de factorizar la matriz A como producto de una matriz triangular inferior por una
triangular superior, a menos de permutaciones. En esta sección profundizamos en este
hecho; concretamente, demostramos el siguiente resultado.
Teorema 2.3.1 (descomposición LU ). Sea A ∈ Mn (R) no singular. Entonces, existen
44
Métodos Numéricos Facultad de Ingenierı́a – IMERL
0 0 1 0
Entonces, siguiendo con la Observación 2.3.1, dada una matriz A ∈ Mn (R), para calcular
B = P A o C = AP en Octave, simplemente podemos escribir
>> B = A(p,:);
>> C = A(:,q);
esto es,
6 −6 6
M A = 0 −1 4 .
0 2 −3
4
45
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Bosquejo de demostración del Teorema 2.3.1. Dada la matriz A ∈ Mn (R), escribimos las
operaciones elementales asociadas al algoritmo de escalerización gaussiana en el paso k-
ésimo:
permutar la fila k-ésima con la fila j-ésima (j > k) se corresponde con multiplicar
a izquierda por una matriz de permutación Pk , que es la matriz identidad con las
filas k y j intercambiadas2 ;
reemplazar una fila dada por una combinación lineal se corresponde con multiplicar
a izquierda por una matriz de multiplicadores Mk .
Además, sabemos que al final del algoritmo de escalerización gaussiana se llega a una
matriz triangular superior, a la que llamaremos U . Ası́, sabemos que existen P1 , . . . , Pn−1
y M1 , . . . , Mn−1 tales que
Mn−1 Pn−1 Mn−2 Pn−2 . . . M2 P2 M1 P1 A = U.
Usando el Ejercicio 2.2 (ver abajo) sucesivamente, la identidad de arriba se puede reescribir
como
Mn−1 M̃n−2 . . . M̃1 Pn−1 Pn−2 . . . P1 A = U,
donde M̃1 , . . . , M̃n−2 son matrices de multiplicadores. Finalmente, definimos
−1
P := Pn−1 Pn−2 . . . P1 , L := Mn−1 M̃n−2 . . . M̃1 .
Tenemos que P es una matriz de permutación, que L es triangular inferior con unos en
la diagonal (ver el Ejercicio 2.3 abajo) y se cumple la identidad deseada,
L−1 P A = U ⇒ P A = LU.
2
Observar que si en el paso k-ésimo no se hace pivoteo, entonces se puede tomar Pk = Id.
46
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ejercicio 2.2. Sea P ∈ Mn (R) de permutación que se obtiene al permutar las filas i y
j de la matriz identidad y sea Mk ∈ Mn (R) de multiplicadores con términos no nulos
debajo de la diagonal en la columna k, con k < i, j. Demostrar que existe M̃ ∈ Mn (R) de
multiplicadores tal que P M = M̃ P . Concretamente, demostrar que M̃ es la matriz que
se obtiene al intercambiar los elementos mik y mjk de la matriz M .
Ejercicio 2.3. Demostrar que el producto de matrices de permutación es de permutación.
Demostrar que el producto de matrices triangulares inferiores con unos en la diagonal es
triangular inferior y con unos en la diagonal. Demostrar que la inversa de una matriz
triangular inferior con unos en la diagonal es triangular inferior y con unos en la diagonal.
Observación 2.3.2 (costo computacional de la factorización
3 LU ). Como realizar el algorit-
mo de escalerización gaussiana tiene un costo de O 2n3 f lops, deducimos que el costo
de computar una descomposición LU de una matriz A es del mismo orden. 4
Observación 2.3.3 (uso de la descomposición LU para resolver un sistema). Supongamos
que queremos resolver un sistema de la forma Ax = b, y ya conocemos una descomposición
de la matriz A como en el Teorema 2.3.1, P A = LU . Multiplicando nuestro sistema a
izquierda por P , obtenemos
LU x = P b.
Llamamos b0 := P b e y := U x. El sistema de arriba se puede pensar como la composición
de dos sistemas triangulares: por un lado el sistema Ly = b0 , y por el otro el sistema U x =
y. Uno puede primero hallar y, y luego utilizar y para hallar x. Como los dos sistemas
involucrados son triangulares, se los puede resolver usando sustitución hacia adelante y
sustitución hacia atrás, respectivamente. El costo de cada una de las resoluciones es de
O(n2 ) f lops, por lo que contar con una descomposición LU de la matriz A nos permite
resolver el sistema Ax = b con O(n2 ) f lops. 4
Observación 2.3.4 (resolver muchos sistemas con una misma matriz). En algunas aplica-
ciones, podemos encontrarnos con la necesidad de resolver muchos sistemas en los que
aparece la misma matriz A ∈ Mn (R), esto es, para vectores b(1) , . . . , b(m) ∈ Rn , tenemos
que hallar x(1) , . . . , x(m) ∈ Rn tales que
47
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Observemos que, si tratamos a una matriz A ∈ Mn (R) como si fuera una matriz llena,
entonces requerimos almacenar n2 números de punto flotante. En cambio, si la tratamos
como una matriz dispersa y tiene m elementos no nulos, entonces requerimos almacenar
3m números de punto flotante. Deducimos, pues, que si 3m n2 , se logra una interesante
reducción de espacio de almacenamiento al tratar a A como dispersa. En Octave, para
indicar que se trate a una matriz como dispersa, se usa la función sparse.
Observación 2.4.1 (inversa de una matriz dispersa). En general, que una matriz invertible
A sea dispersa no implica que A−1 también lo sea. 4
Una clase de matrices dispersas muy especial y que aparece en múltiples aplicaciones es
la de aquellas matrices cuyos elementos no nulos se concentran alrededor de la diagonal.
3
A veces, se utiliza el inexistente término esparsa para nombrar a estas matrices; esta palabra parece
derivar del inglés sparse.
48
Métodos Numéricos Facultad de Ingenierı́a – IMERL
a las matrices que tienen ancho de banda igual a 1 se las llama tridiagonales.
>> x = A \ b;
49
Métodos Numéricos Facultad de Ingenierı́a – IMERL
1. kuk ≥ 0 ∀u ∈ Rn , y kuk = 0 ⇔ u = 0;
2. kλ uk = |λ| kuk ∀λ ∈ R, ∀u ∈ Rn ;
50
Métodos Numéricos Facultad de Ingenierı́a – IMERL
De cursos anteriores de álgebra lineal, sabemos que el espacio de matrices Mm×n (R), con
las operaciones habituales de suma de matrices y de producto de matrices por números
reales, tiene estructura de espacio vectorial. Para definir normas, puede ser tentador iden-
tificar matrices en Mm×n (R) con vectores en Rmn y utilizar las normas `p o `∞ de vectores
que discutimos arriba4 . Aquı́ nos concentramos en normas de matrices definidas de otra
forma, las llamadas normas de operador, que están asociadas a identificar las matrices en
Mm×n (R) con transformaciones lineales Rn → Rm .
Definición 2.5.1 (norma inducida o norma operador). Sea k · kv una norma vectorial.
Definimos la norma matricial inducida por k · kv (o norma operador asociada a
k · kv ), k · k : Mm×n (R) → R, mediante
kAxkv
kAk := máx . (2.5)
x6=0 kxkv
Ejercicio 2.5. a) Verificar que las normas inducidas son efectivamente normas. Esto es,
demostrar que la norma dada por (2.5) verifica las tres condiciones que mencionamos
al comienzo de esta sección.
kAxkv
kAk = máx = máx kAxkv = máx kAxkv .
x6=0 kxkv kxkv ≤1 kxkv =1
Esto muestra que el valor de la norma inducida de una matriz A por una cierta norma
vectorial k · kv es igual al máximo que logra “estirar” (según la norma k · kv ) la matriz
A a los vectores de norma k · kv igual a 1.
4
El caso p = 2 da lugar a la llamada norma de Frobenius, que es de interés en varias aplicaciones.
51
Métodos Numéricos Facultad de Ingenierı́a – IMERL
La siguiente proposición muestra que la definición anterior aplica a las normas matriciales
inducidas.
Proposición 2.5.1 (compatibilidad y submultiplicatividad de las normas operador). Sea
k · k la norma matricial inducida por una norma vectorial k · kv . Se cumplen las desigual-
dades
kAxkv ≤ kAk kxkv ∀A ∈ Mm×n (R), ∀x ∈ Rn ,
y
kABk ≤ kAkkBk ∀A ∈ Mm×n (R), B ∈ Mn×p (R). (2.6)
52
Métodos Numéricos Facultad de Ingenierı́a – IMERL
A partir de este punto, cuando escribamos normas no haremos referencia explı́cita a qué
tipo de norma estamos considerando (matricial o vectorial), sino que esto se deduce del
argumento que tome la norma en cada caso. Asimismo, cuando escribamos una norma
matricial, por defecto asumimos que es una norma operador asociada a una cierta norma
vectorial.
Esto es, la norma matricial inducida por la norma `∞ es igual al máximo de las normas
`1 de las filas de la matriz.
Demostración. Sea A = (aij ) ∈ Mm×n (R). Si A es la matriz nula, entonces (2.7) vale
trivialmente. Asumamos A no es nula y tomemos un vector x ∈ Rn con kxk∞ = 1, lo que
implica que
|xj | ≤ máx |xk | = kxk∞ = 1 ∀j = 1, . . . , n.
k
Usamos esta propiedad junto a la desigualdad triangular, y la definición de norma `∞
vectorial para obtener
n
X n
X n
X
kAxk∞ = máx |(Ax)i | = máx aij xj ≤ máx |aij xj | ≤ máx |aij |.
i=1,...,m i=1,...,m i=1,...,m i=1,...,m
j=1 j=1 j=1
Por lo tanto, utilizando la segunda parte del Ejercicio 2.5, llegamos a la desigualdad
n
X
kAk∞ = máx kAxk∞ ≤ máx |aij |. (2.8)
kxk∞ =1 i=1,...,m
j=1
Para probar que se cumple (2.7), alcanza con probar que la cota hallada se alcanza. Es
decir, basta con hallar un vector y ∈ Rn con kyk∞ = 1 tal que la desigualdad de arriba
sea de hecho una igualdad. Para este fin, consideremos un ı́ndice i0 ∈ 1, . . . , m tal que la
fila i0 de la matriz A tiene la mayor norma `1 entre todas las filas de A,
n
X n
X
|ai0 j | = máx |aij |,
i=1,...,m
j=1 j=1
53
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Sea A ∈ Mm×n (R). Ası́ como la norma matricial de A inducida por la norma `∞ es igual al
máximo de las normas `1 de las filas de A, la norma matricial de A inducida por la norma
`1 es igual al máximo de las normas `1 de las columnas de A. Dejamos la demostración
de esta propiedad como ejercicio.
Ejercicio 2.6. Sea k·k1 la norma matricial inducida por la norma `1 vectorial. Demostrar
que vale la igualdad
m
X
kAk1 = máx |aij | ∀A = (aij ) ∈ Mm×n (R).
j=1,...,n
i=1
La norma matricial inducida por la norma euclı́dea también tiene una caracterización
importante. En ella aparecen los llamados valores singulares de la matriz. Vamos a pro-
fundizar en ellos cuando tratemos problemas de mı́nimos cuadrados.
Proposición 2.5.4 (norma inducida por la norma euclı́dea). Sea k·k2 la norma matricial
inducida por la norma `2 vectorial. Para toda A ∈ Mm×n (R), se cumple que
p
kAk2 = λ1 , (2.9)
√
donde λ1 es el mayor valor propio de At A. El número σ1 (A) := λ1 se llama el primer
valor singular de la matriz A.
Observación 2.5.3 (valores propios de At A). Antes de demostrar la proposición, observe-
mos que el enunciado tiene sentido: que At A tiene valores propios reales y son no negativos.
Dada A ∈ Mm×n (R) es sencillo verificar que la matriz At A ∈ Mn (R) es simétrica; por
lo tanto, el Teorema Espectral para matrices reales y simétricas nos asegura que At A es
diagonalizable. Además, At A es semidefinida positiva, pues para todo x ∈ Rn tenemos
(At Ax) · x = xt At Ax = (Ax)t (Ax) = kAxk22 ≥ 0. (2.10)
Esto quiere decir que para toda A ∈ Mm×n (R), todos los valores propios de At A son
mayores o iguales que 0. Por lo tanto, tiene sentido tomar raı́z cuadrada en (2.9). 4
54
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por lo tanto, usando propiedades elementales del producto de matriz por vector, la defi-
nición de los coeficientes αi y el hecho de que λ1 ≥ λi para todo i, tenemos
n
X n
X n
X
t t t t t t t
x A Ax = x (A Ax) = x αi A Avi = x αi λi vi = αi2 λi ≤ λ1 .
i=1 i=1 i=1
Tomando raı́z cuadrada en la desigualdad anterior√ y tomando máximo entre todos los
vectores con kxk2 = 1, deducimos que kAk2 ≤ λ1 .
√
Para probar
√ que vale kAk 2 ≥ λ1 , basta con encontrar un vector y tal que kyk2 = 1 y
kAyk2 = λ1 . Elegimos y como un vector propio de At A asociado al mayor valor propio
λ1 . Usando la definición de norma euclı́dea, tenemos que
p p p p p p
kAyk2 = yt At Ay = yt λ1 y = λ1 yt y = λ1 kyk2 = λ1 .
55
Métodos Numéricos Facultad de Ingenierı́a – IMERL
56
Métodos Numéricos Facultad de Ingenierı́a – IMERL
57
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 2.1: Izquierda: solución computada x̄ y solución exacta x del sistema (2.12). De-
recha: un sistema que, al resolverlo y obtener residuo pequeño, podemos asegurar que el
error es pequeño.
58
Métodos Numéricos Facultad de Ingenierı́a – IMERL
kAxkv
Demostración. Por el Ejercicio 2.5, ya sabemos que kAk = máxx6=0 kxkv
, por lo que resta
demostrar que
1
kA−1 k = .
mı́nx6=0 kAxk
kxkv
v
59
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por otra parte, usando forma análoga las identidades equivalentes b = Ax y x = A−1 b,
obtenemos
kbk
≤ kxk ≤ kA−1 k kbk (2.15)
kAk
Combinando (2.14) y (2.15), deducimos
1 krk kek krk
−1
≤ ≤ kAk kA−1 k
kAk kA k kbk kxk kbk
y usando la Definición 2.5.3 concluimos
1 krk kek krk
≤ ≤ κ(A) . (2.16)
κ(A) kbk kxk kbk
Notemos que el término en el centro arriba es un error relativo en x, mientras que el
krk
término kbk corresponde a un residuo relativo. La identidad (2.16) nos dice entonces que
este residuo relativo y el error relativo al aproximar la solución de nuestro sistema Ax = b
están relacionados mediante el número de condición de la matriz A. En otras palabras,
si κ(A) no es muy grande, residuos relativos pequeños se corresponden con
errores relativos pequeños; en cambio, si κ(A) es grande, el hecho de que el
residuo relativo sea pequeño no nos permite asegurar que el error relativo lo
sea.
60
Métodos Numéricos Facultad de Ingenierı́a – IMERL
kδx k kδb k
≤ κ(A) .
kxk kbk
61
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(A + E) x̄ = b, (2.18)
donde E ∈ Mn (R) es una matriz cuyos elementos son (en valor absoluto) del orden de
los errores de redondeo al representar la matriz A.
Analicemos el enunciado del Teorema 2.5.6. El hecho de que los elementos de E sean del
orden de los errores de redondeo al representar la matriz A significa que, para cualquier
norma matricial, tenemos
kEk ≤ ρ εM kAk
para alguna constante ρ > 0, que podemos esperar cumpla ρ . 10. Por lo tanto, ese
teorema nos garantiza estar en una situación en la que nuestra solución computacional x̄
es la solución exacta a un sistema al que le hemos perturbado la matriz A, la perturbación
es δA = E, y tenemos un cierto control en norma sobre esta perturbación. Notamos que
δx = x̄ − x = e y que por lo tanto la fórmula (2.17) se puede reescribir como
kek kEk
≤ κ(A) ≤ κ(A) ρ εM .
kx̄k kAk
krk kEk
≤ ≤ ρ εM .
kAk kx̄k kAk
7
Para entender por qué esta es una buena normalización del residuo, puede servir pensar en el sistema
como si A, x y b tuviesen dimensiones (de largo, de tiempo, etc). Las dimensiones del residuo tienen que
ser iguales que las del vector b, que a su vez deben ser iguales al producto de las dimensiones de A por
las de x.
62
Capı́tulo 3
Interpolación polinomial
0
y
−1
−2
−3
−4
0 0,5 1 1,5 2 2,5 3
x
Figura 3.1: Cualquiera de las curvas dibujadas pasa por los tres puntos marcados.
63
Métodos Numéricos Facultad de Ingenierı́a – IMERL
En forma similar, dada una tabla de evaluaciones de una cierta función (desconocida)
f , nos puede interesar aproximar la derivada o la integral de esta función f . Para
ello, podemos nuevamente construir una interpolante de los datos de la tabla y
utilizarla para derivarla o integrarla analı́ticamente.
64
Métodos Numéricos Facultad de Ingenierı́a – IMERL
para todo i = 0, . . . , n. La primer pregunta que surge es qué requerimientos tiene sentido
pedirle a p: ¿puede ser un polinomio de cualquier grado?
Razonemos intuitivamente. Si tenemos un único punto (x0 , y0 ), entonces hay un único
polinomio de grado 0 que cumple lo estipulado, la función constante p(x) = y0 , y ob-
viamente hay infinitos polinomios de grado mayor o igual a 1 cuyos gráficos pasan por
(x0 , y0 ). Si tenemos dos puntos {(x0 , y0 ), (x1 , y1 )}, entonces en general –salvo que y0 = y1 –
no podemos esperar que el gráfico de una función constante pase por esos dos puntos. En
cambio, sı́ sabemos que existe una única recta por los puntos (x0 , y0 ), (x1 , y1 ) y podemos
considerar la función lineal cuyo gráfico sea dicha recta: para n + 1 = 2 puntos, hay una
única función lineal que cumple lo estipulado. Si consideramos polinomios de grado más
alto, entonces tenemos demasiados grados de libertad y no tenemos una solución única.
Del razonamiento de arriba parece deducirse que, si esperamos que nuestro problema de
interpolación por n + 1 puntos tenga solución y que ésta sea única, entonces debemos
requerir que el polinomio p sea de grado n. El siguiente teorema demuestra que esta
intuición es correcta.
Teorema 3.2.1 (polinomio interpolante). Sea n ≥ 0 un número entero. Dados n + 1
puntos en el plano {(xi , yi )}i=0,...,n con xi 6= xj si i 6= j, existe un único polinomio pn de
grado menor o igual que n tal que pn (xi ) = yi para todo i = 0, . . . , n.
Definición 3.2.1 (polinomio interpolante). Al polinomio definido por el teorema anterior
lo llamamos el polinomio interpolante por los puntos {(xi , yi )}i=0,...,n . 4
65
Métodos Numéricos Facultad de Ingenierı́a – IMERL
y an ∈ R es una incógnita a determinar. Para hallar an , nos resta imponer que pn (xn ) = yn :
evaluando (3.1) en xn y usando la expresión para q, tenemos
Ası́, tenemos determinado un polinomio pn de grado menor o igual a n que cumple pn (xi ) =
yi para todo i = 0, . . . , n. La unicidad de pn se deduce de que todo polinomio de grado
menor o igual a n se puede escribir como (3.1) y del hecho de que el an que hallamos es
único.
El teorema anterior nos da condiciones para que el problema de interpolar datos con una
función polinomial tenga solución única. A continuación, exploramos distintas formas de
hallar este polinomio interpolante.
Remarcamos que aquı́ las incógnitas son los coeficientes c0 , . . . , cn del polinomio pn ; los
puntos (x0 , y0 ), . . . , (xn , yn ) son datos de nuestro problema. En forma matricial, el sistema
anterior se puede escribir como
1 x0 x20 . . . xn0 c0 y0
1 x1 x2 . . . xn c1 y1
1 1
.. .. = .. . (3.5)
.. .. ..
. . . . . .
1 xn x2n . . . xnn cn yn
66
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Definición 3.2.2 (matriz de Vandermonde). A toda matriz que presenta una progresión
geométrica en cada fila, tal como la matriz que aparece en el sistema (3.5), se le llama
matriz de Vandermonde. 4
Observación 3.2.1 (invertibilidad de matrices de Vandermonde). Si xi 6= xj cuando i 6= j,
entonces la matriz de Vandermonde asociada es invertible. Esto es una consecuencia inme-
diata del Teorema 3.2.1, que nos asegura la existencia y unicidad del polinomio interpolan-
te. Si escribimos el polinomio interpolante como en (3.4), entonces el Teorema 3.2.1 puede
interpretarse como que el sistema (3.5) sea compatible determinado. En consecuencia, la
matriz de dicho sistema, que es la matriz de Vandermonde, es invertible. 4
Ejemplo 3.2.1 (forma de Vandermonde). Tenemos la siguiente tabla con tres datos:
i xi yi
0 1/4 -3/4
1 1 -1
2 5/2 3/2
Sabemos que existe un único polinomio p2 de grado menor o igual a 2 tal que su gráfico
pasa por esos tres puntos. Si escribimos p2 como en (3.4), nuestras incógnitas son los
coeficientes c0 , c1 , c2 , y las podemos hallar resolviendo el sistema (3.5), que en este ejemplo
toma la forma
1 1/4 1/16 c0 −3/4
1 1 1 c1 = −1 .
1 5/2 25/4 c2 3/2
Podemos resolver este sistema en forma sencilla usando cualquiera
de las técnicas que
c0 −4/9
discutimos el capı́tulo anterior, y llegamos a la solución c1 = −13/9. Deducimos
c2 8/9
entonces que el polinomio interpolante por los datos de la tabla es
4 13 8
p2 (x) = − − x + x2 . (3.6)
9 9 9
Se puede verificar que este polinomio verifica p2 (xi ) = yi para i = 0, 1, 2, y su gráfico es
el que está representado por la curva azul en la Figura 3.1. 4
67
Métodos Numéricos Facultad de Ingenierı́a – IMERL
>> x = linspace(0,1,21);
>> V = vander(x);
>> condest(V)
ans = 6.7512e+16
En la primera lı́nea, usamos el comando linspace para definir un vector fila x que tiene
21 elementos equiespaciados en el intervalo [0, 1], esto es, x= [0, 1/20, . . . , 1]. La segunda
lı́nea usa la función de Octave vander para definir la matriz de Vandermonde V asociada
a los nodos x(i), y la tercera estima el número de condición de V asociado a la norma
`1 (recordar la Definición 2.5.3). Observamos que el resultado es bastante desalentador:
nuestra estimación es que κ(V) es del orden de 1016 . Esto implica, por ejemplo y teniendo
en cuenta lo que hemos discutido en la Sección 2.5.4, que si quisiéramos interpolar por
21 puntos cuyas abscisas sean las del vector x y usáramos la forma de Vandermonde
resolviendo el sistema (3.5) mediante escalerización gaussiana con pivoteo parcial, nuestra
mejor cota para el error relativo en los coeficientes serı́a del orden de 1. En otras palabras,
¡podrı́amos obtener un error del 100 % en los coeficientes del polinomio!
Notar, sin embargo, que esto no implica que nuestra evaluación del polinomio interpolante
sea mala, ya que al hacer escalerización gaussiana con pivoteo parcial el residuo es del
orden de εM independientemente del condicionamiento de la matriz. 4
Observación 3.2.3 (agregar datos en la forma de Vandermonde). Otra limitación de la
forma de Vandermonde surge si queremos agregar puntos por los que interpolar.
Supongamos que tenemos la misma tabla de datos que en el Ejemplo 3.2.1, hallamos p2
como allı́ se indica, pero ahora disponemos de un nuevo dato adicional (x3 , y3 ). Esto es,
queremos hallar un polinomio p3 de grado menor o igual a 3 que interpole por los puntos
(x0 , y0 ), . . . , (x3 , y3 ). Si usamos la forma de Vandermonde, entonces el trabajo que hicimos
para hallar p2 no servirá de mucho, porque vamos a tener que rearmar el sistema
(3.5) pero ahora con la matriz de Vandermonde asociada a los 4 puntos x0 , x1 , x2 , x3 . En
general, al usar (3.4) no existe una relación directa entre la solución del sistema que nos
da los coeficientes de p2 con el que nos dé los coeficientes de p3 , por lo que tendremos que
empezar de cero nuevamente. 4
68
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(x0 , 0), . . . , (xk−1 , 0), (xk , 1), (xk+1 , 0), . . . , (xn , 0).
En efecto, basta con observar que el polinomio dado en (3.7) verifica la Definición 3.2.3 y
usar la unicidad del polinomio interpolante dada por el Teorema 3.2.1. 4
Hallar la funciones de base de Lagrange puede ser laborioso, pero una vez calculadas
{Lkn }k=0,...,n , el polinomio interpolante se puede expresar en forma sumamente sencilla.
Por lo tanto, el lado derecho de (3.8) coincide con el polinomio interpolante por
{(xi , yi )}i=0,...,n .
Ejemplo 3.2.2 (forma de Lagrange). Volvemos al Ejemplo 3.2.1 y a la tabla de tres datos
allı́ presentada. Usando la expresión (3.7), encontramos que la base de Lagrange asociada
69
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(x − 1)(x − 5/2) 16
L02 (x) = = (x − 1)(x − 5/2),
(1/4 − 1)(1/4 − 5/2) 27
(x − 1/4)(x − 5/2) 8
L12 (x) = = − (x − 1/4)(x − 5/2),
(1 − 1/4)(1 − 5/2) 9
(x − 1/4)(x − 1) 8
L22 (x) = = (x − 1/4)(x − 1).
(5/2 − 1)(5/2 − 1/4) 27
La Figura 3.2 muestra estas tres funciones de base. Por lo tanto, usando los valores de las
coordenadas y de la tabla de datos, deducimos que el polinomio interpolante es
2
X 3 3
p2 (x) = yk Lk2 (x) = − L02 (x) − L12 (x) + L22 (x).
k=0
4 2
Dejamos para el lector verificar que esta expresión coincide con (3.6). 4
1,5
0,5
0
y
−0,5
−1 L02
L12
−1,5 L22
−2
0 0,5 1 1,5 2 2,5 3
x
70
Métodos Numéricos Facultad de Ingenierı́a – IMERL
entonces la evaluación de (3.7) para hallar la base correspondiente puede ser costosa.
Además, si queremos evaluar pn en algún punto que no sea x0 , . . . , xn , entonces será más
costoso computacionalmente hacerlo si lo tenemos escrito en la base de Lagrange que si
lo tenemos escrito en la base monomial. En forma análoga, la forma de Lagrange puede
ser más difı́cil de manipular a la hora de integrar o derivar el polinomio interpolante, lo
que es un uso relevante de esta herramienta (ver la Sección 3.1).
Por otra parte, al igual que ocurre con la forma de Vandermonde, la forma de Lagrange no
es demasiado amigable con agregar nuevos puntos para interpolar: si se agrega un nuevo
punto, será necesario recalcular la nueva base. 4
Ejercicio 3.1. Sean los coeficientes {ai }i=0,...,n dados por (3.3). Verificar que el lado
derecho de (3.9) es efectivamente una expresión para el polinomio interpolante por los
puntos {(xi , yi )}i=0,...,n .
[Sugerencia: evaluar el lado derecho en x0 , . . . , xn y usar la unicidad del polinomio interpolante.]
Observación 3.2.7 (sobre la demostración del Teorema 3.2.1). Nuestra demostración del
Teorema 3.2.1 está relacionada con la forma de Newton y motivada por ella. Ese teorema se
puede demostrar de varias formas alternativas. En particular, referimos a [Hea02, Sección
7.3.1] para una demostración basada en la forma de Vandermonde y a [SM03, Teorema
6.1] para una inspirada en la forma de Lagrange. 4
71
Métodos Numéricos Facultad de Ingenierı́a – IMERL
El Algoritmo 3.1 muestra cómo esta expresión puede ser implementada de forma muy
sencilla.
Algoritmo 3.1: Pseudo-código: algoritmo de Horner
Datos: a0 , . . . , an ∈ R, x ∈ R
Resultado: y = pn (x)
y ← an ;
para j = n − 1 : −1 : 0 hacer
y ← aj + (x − xj )y;
fin
4
Observación 3.2.9 (diferencias divididas). Una vez que tenemos determinados los coefi-
cientes {ai }i=0,...,n , el algoritmo de Horner nos permite evaluar pn de forma eficiente. Para
determinar estos coeficientes, bien podrı́amos utilizar (3.3) directamente, aunque en esta
expresión también hay una “recursión escondida”.
Esta recursión se explota en el llamado método de diferencias divididas, que descri-
bimos a continuación. Se parte2 de f [xi ] = f (xi ) = yi y se toma, para i = 1, . . . , n,
f [xk+1 , xk+2 , . . . , xk+i ] − f [xk , xk+1 , . . . , xk+i−1 ]
f [xk , xk+1 , . . . , xk+i ] := , k = 0, . . . , n − i.
xk+i − xk
Una vez halladas estas diferencias divididas, simplemente se tiene
a0 = f [x0 ], a1 = f [x0 , x1 ], . . . , an = f [x0 , x1 , . . . , xn ].
4
Ejemplo 3.2.3 (forma de Newton). Volvamos al Ejemplo 3.2.1/3.2.2, y expresemos el
polinomio interpolante p2 en la forma de Newton.
En primer lugar, hallamos los coeficientes a0 , a1 , a2 mediante el método de diferencias
divididas. Recordando los valores de la tabla de datos de este ejemplo, tenemos
f [x0 ] = −3/4, f [x1 ] = −1, f [x2 ] = 3/2 (i = 0, k = 0),
f [x1 ] − f [x0 ] −1 − (−3/4) 1
f [x0 , x1 ] = = =− (i = 1, k = 0),
x1 − x0 1 − 1/4 3
f [x2 ] − f [x1 ] 3/2 − (−1) 5
f [x1 , x2 ] = = = (i = 1, k = 1),
x2 − x1 5/2 − 1 3
f [x1 , x2 ] − f [x0 , x1 ] 5/3 − (−1/3) 8
f [x0 , x1 , x2 ] = = = , (i = 2, k = 0).
x2 − x0 5/2 − 1/4 9
2
A no asustarse con la expresión f [·]: ¡es solamente notación!
72
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Queda como ejercicio para el lector verificar que esta expresión coincide con (3.6). 4
Definición 3.3.1 (norma del supremo). Sea f : D → R una función. Definimos la norma
del supremo de f en D, que escribimos como kf kL∞ (D) , mediante
73
Métodos Numéricos Facultad de Ingenierı́a – IMERL
en := f − pn .
Notar que en : [a, b] → R es una función: para cada punto x ∈ [a, b] tenemos un valor
en (x). Teniendo en cuenta la Definición 3.3.1, nos interesa estimar ken kL∞ ([a,b]) .
f (n+1) (γx )
en (x) = f (x) − pn (x) = (x − x0 )(x − x1 ) . . . (x − xn ). (3.11)
(n + 1)!
La fórmula (3.11) nos dice que el error de interpolación en un punto dado x ∈ [a, b]
depende de la derivada (n + 1)-ésima de f en un cierto punto (desconocido y que depende
de x) del intervalo [a, b] y del producto de las distancias (con signo) entre x y los nodos
de interpolación x0 , . . . , xn . Si nos interesa acotar el error de interpolación polinomial,
entonces debemos estimar la magnitud de estos términos: volvemos a este punto en la
Observación 3.3.2 más adelante. Por el momento, y para alivianar un poco la notación,
introducimos la siguiente definición.
Definición 3.3.2 (polinomio nodal). Sean los puntos x0 < x1 < . . . < xn en el intervalo
[a, b]. Definimos el polinomio nodal por x0 , . . . , xn como ωn : [a, b] → R,
n
Y
ωn (x) := (x − xi ).
i=0
Demostración del Teorema 3.3.1. Fijemos un punto x ∈ [a, b]. Es claro que si x = xi para
algún i entonces f (xi ) − pn (xi ) = 0 y vale (3.11), por lo que podemos suponer que x 6= xi
para todo i = 0, . . . , n.
5
Observar que aquı́ estamos definiendo el error con el signo opuesto a como lo hicimos en el Capı́tulo
1. De todas formas, nos interesa la magnitud (valor absoluto) de este error y no su signo.
74
Métodos Numéricos Facultad de Ingenierı́a – IMERL
ωn (t)
G(t) := en (t) − en (x) . (3.12)
ωn (x)
(n + 1)!
G(n+1) (t) = f (n+1) (t) − en (x) ,
ωn (x)
y evaluando en γx ,
(n + 1)!
0 = G(n+1) (γx ) = f (n+1) (γx ) − en (x) .
ωn (x)
La utilidad principal del Teorema 3.3.1 radica en poder acotar el error de interpolación.
Corolario 3.3.2 (estimación del error de interpolación polinomial). Bajo las mismas
hipótesis del Teorema 3.3.1, tenemos las estimaciones
75
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Demostración. Ambas estimaciones son consecuencia directa de (3.11). Dado x ∈ [a, b],
tomando valor absoluto en (3.11) y acotando |f (n+1) (γx )| ≤ kf (n+1) kL∞ ([a,b]) , deducimos
(3.13). Luego, como además vale |ωn (x)| ≤ kωn kL∞ ([a,b]) , tomando supremo en x ∈ [a, b]
obtenemos (3.14).
f 0 (x) = cos2 (x) − sen2 (x), f 00 (x) = −4 sen(x) cos(x), f 000 (x) = −4(cos2 (x) − sen2 (x)), . . . .
Acotando (groseramente) el valor absoluto de las funciones seno y coseno por 1, deducimos
que para todo x ∈ [0, 1] valen
76
Métodos Numéricos Facultad de Ingenierı́a – IMERL
lo que podemos escribir como kf (n) kL∞ ([0,1]) ≤ 2n ∀n ≥ 1. Por otra parte, como esta-
mos trabajando en el intervalo [0, 1], para todo x ∈ [0, 1] podemos acotar (nuevamente,
groseramente)
|x − xi | ≤ 1 ∀i = 0, . . . , n,
lo que implica
n
Y
kωn kL∞ ([0,1]) = máx |ωn (x)| = máx |x − xi | ≤ 1.
x∈[0,1] x∈[0,1]
i=0
Si estamos en una situación como la del Ejemplo 3.3.1, en que podemos estimar
kωn kL∞ ([0,1]) de forma razonable, y las normas del supremo de las derivadas de
la función f no crecen “demasiado rápido” respecto a n, entonces es de esperar
que al aumentar la cantidad de puntos por los que interpolamos, el error de
interpolación decrezca. En estos casos, vale la pena hacer una interpolación polinomial
de grado alto. 4
77
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 3.4: Fenómeno de Runge: interpolar la función de Runge sobre puntos equiespa-
ciados.
¿A qué se debe tan mal comportamiento del polinomio interpolante? Recordemos los tres
puntos que marcamos en la Observación 3.3.2. Se puede demostrar –es una cuenta larga y
78
Métodos Numéricos Facultad de Ingenierı́a – IMERL
1. Mantenerse con interpolaciones de grado bajo. Esto logra evitar las oscila-
ciones que observamos en la Figura 3.4. Como una interpolación de grado bajo en
un intervalo “grande” es poco precisa, la estrategia habitual consiste en dividir el
intervalo en subintervalos pequeños y en cada uno de estos subintervalos aplicar
una interpolación de grado bajo. Vamos a explorar en detalle esta estrategia en la
Sección 3.4.
2. Elegir los nodos de interpolación de forma más astuta, de modo de lograr que
kωn kL∞ ([−1,1]) sea lo menor posible. Comentamos muy brevemente esta alternati-
va. Los nodos de Chebyschev en el intervalo [−1, 1] son
πi
xi = − cos , i = 0, . . . , n,
n
y tienen la propiedad de minimizar
n
Y
máx |x − xi | = kωn kL∞ ([−1,1]) .
x∈[−1,1]
i=0
La Figura 3.5 nos muestra las interpolaciones polinomiales p5 , p10 , p20 de la fun-
ción de Runge utilizando los nodos de Chebyschev. Notemos que no se observan las
oscilaciones presentes en la Figura 3.4 y que efectivamente una buena elección de
los nodos por los que interpolamos nos puede permitir estabilizar estas oscilacio-
nes. De hecho, se puede demostrar que con esta elección se tiene la convergencia
ken kL∞ ([−1,1]) → 0 con n → ∞. Referimos a [QS06, Sección 3.1.2] para más detalles.
79
Métodos Numéricos Facultad de Ingenierı́a – IMERL
p [xi−1 ,xi ]
(xi ) = p [xi ,xi+1 ]
(xi ) y p0 [xi−1 ,xi ]
(xi ) = p0 [xi ,xi+1 ]
(xi ). (3.15)
Por lo tanto, si queremos que el polinomio interpolante sea derivable, debemos prescribir
de alguna forma el valor de p0 (xi ); llamemos di a ese valor, que por el momento asumimos
conocido. Si nos restringimos a un intervalo [xi , xi+1 ] (i = 0, . . . , n − 1), requerimos que p
80
Métodos Numéricos Facultad de Ingenierı́a – IMERL
cumpla 4 condiciones:
Como tenemos 4 restricciones, si queremos tener un problema con solución única parece
natural que lo podamos lograr con 4 grados de libertad. Por lo tanto, en cada intervalo
[xi , xi+1 ] (i = 0, . . . , n − 1) el polinomio p debe ser un polinomio cúbico.
Más en general, al hacer interpolaciones polinomiales a trozos, las funciones
interpolantes se suelen tomar como polinomios de grado impar en cada subin-
tervalo. En esta sección nos restringimos a funciones lineales y cúbicas a trozos.
81
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ejercicio 3.2. El Algoritmo 3.2 computa todas las diferencias divididas, lo que es innece-
sario si solamente nos interesa evaluar L en un punto. Sin embargo, esto puede ser útil si
queremos evaluar L en varios puntos a la vez. Modificando el Algoritmo 3.2, implementar
un código que permita evaluar L en un vector de puntos v ∈ Rm de forma eficiente, asu-
miendo que las entradas de v están ordenadas de forma creciente y pertenecen al intervalo
[x0 , xn ].
En forma similar a la Sección 3.3, cabe preguntarse qué tan bien podemos aproximar a
una cierta función f si usamos su interpolación lineal a trozos. Concretamente, si tenemos
que yi = f (xi ) para todo i = 0, . . . , n, y construimos la interpolante lineal a trozos L,
queremos estimar la discrepancia entre f y L. Nuestra herramienta principal es el Teorema
3.3.1.
En efecto, supongamos que f es de clase C 2 , tomemos n + 1 puntos x0 < x1 < . . . < xn , y
L la interpolante lineal a trozos por {(xi , f (xi )}i=0,...,n . Consideremos un punto x ∈ [x0 , xn ]
fijo: si x = xi para algún i = 0, . . . , n, entonces el error de interpolación es nulo. Si x no
coincide con ninguno de los nodos, entonces x ∈ (xi , xi+1 ) para algún i = 0, . . . , n − 1.
Observemos que en ese intervalo estamos haciendo una interpolación lineal (n = 1) de f ,
por lo que podemos usar (3.11): existe un γx ∈ (xi , xi+1 ) tal que
f 00 (γx )
f (x) − L(x) = (x − xi )(x − xi+1 ). (3.17)
2
Esta fórmula nos da una representación del error de interpolación; si queremos estimar
su magnitud, entonces debemos acotar el lado derecho de la igualdad de forma análoga a
como hicimos en la demostración del Corolario 3.3.2.
Teorema 3.4.1 (estimación del error de interpolación lineal a trozos). Sean n + 1 puntos
x0 < x1 < . . . < xn , f una función de clase C 2 en el intervalo [x0 , xn ], y L la interpolante
lineal a trozos por los puntos {(xi , f (xi ))}i=0,...,n . Entonces, si x ∈ [xi , xi+1 ] para i ∈
{0, . . . , n − 1}, se tiene
kf 00 kL∞ ([xi ,xi+1 ])
|f (x) − L(x)| ≤ (x − xi ) (xi+1 − x). (3.18)
2
82
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Además,
kf 00 kL∞ ([x0 ,xn ])
kf − LkL∞ ([x0 ,xn ]) ≤ máx |xj+1 − xj |2 . (3.19)
8 j=0,...,n−1
Demostración. Sea x ∈ [xi , xi+1 ] para algún i ∈ {0, . . . , n − 1}. Tomando valor absoluto
en (3.17), y como γx ∈ (xi , xi+1 ) podemos acotar
alcanza su máximo en x = xi +x2 i+1 , el punto medio del intervalo [xi , xi+1 ], y este máximo
2
vale |xi+14−xi | . Por lo tanto, si x ∈ [xi , xi+1 ], tenemos
2i
xi := −1 + , i = 0, . . . , n,
n
y construimos la interpolante lineal a trozos de f por estos nodos. Esto es, en la Definición
3.4.1 tomamos yi = f (xi ) para i = 0, . . . , n. La Figura 3.6 nos muestra el resultado
tomando n + 1 = 21 puntos, y se contrasta notablemente con el panel derecho de la
Figura 3.4: la interpolación lineal a trozos evita que la función interpolante oscile.
También podemos estimar el error de interpolación usando el Teorema 3.4.1. La derivada
segunda de la función de Runge es
50 (75x2 − 1)
f 00 (x) = ,
(1 + 25x2 )3
83
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por otra parte, estamos tomando nodos equiespaciados, por lo que todos los subintervalos
tienen la misma longitud: tenemos el intervalo [−1, 1] partido en n = 20 subintervalos,
por lo que tenemos
1
|xi+1 − xi | = ∀i = 0, . . . , n − 1.
10
Reemplazando en (3.19), obtenemos que, al realizar una interpolación lineal a trozos de
la función de Runge con 21 nodos equiespaciados, el error en norma del supremo se puede
acotar mediante
50 1 1
kf − LkL∞ ([−1,1]) ≤ 2
= .
8 10 16
Más en general, si en lugar de 21 usáramos n + 1 nodos equiespaciados, entonces cada
subintervalo tendrı́a longitud
2
|xi+1 − xi | = ∀i = 0, . . . , n − 1.
n
y llegarı́amos a la cota del error
2
50 2 25
kf − LkL∞ ([−1,1]) ≤ = 2.
8 n n
Notemos que kf − LkL∞ ([−1,1]) → 0 cuando n → ∞: a diferencia de lo que ocurrı́a en la
Figura 3.4, ahora tenemos que la sucesión de interpolantes lineales a trozos por n puntos
equiespaciados converge a la función. 4
84
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Es conveniente tratar este problema en cada subintervalo por separado. Para probar que
existe un único polinomio cúbico en el intervalo [xi , xi+1 ] que cumple
3hi s2 − 2s3
3
hi − 3hi s2 + 2s3
p(x) = yi+1 + yi
h3i h3i
(3.21)
s2 (s − hi ) s(s − hi )2
+ di+1 + di .
h2i h2i
s := x − xi .
Ejercicio 3.3. Verificar que el polinomio dado por (3.21) cumple las condiciones (3.20).
La diferencia entre las tres alternativas que analizamos a continuación radica en cómo
tomamos las pendientes d0 , . . . , dn . Dependiendo del problema que nos interese, éstas
podrı́an ser conocidas o no. En caso de que no lo sean, tenemos la libertad de elegirlas a
nuestro criterio.
7
Notar que, como el sistema resultante es 4 × 4, no tenemos los inconvenientes que comentamos en las
observaciones 3.2.2 y 3.2.3.
85
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Interpolación de Hermite
yi = f (xi ), di = f 0 (xi ), ∀i = 0, . . . , n.
Como las pendientes d0 , . . . , dn son conocidas, simplemente debemos reemplazar los va-
lores correspondientes en (3.21) y logramos construir una interpolante cúbica a trozos
p : [x0 , xn ] → R, que cumple p(xi ) = f (xi ), p0 (xi ) = f 0 (xi ) para todo i = 0, . . . , n.
Definición 3.4.2 (interpolante cúbica a trozos de Hermite). A la función p : [x0 , xn ] → R
construida como se describe arriba se le llama la interpolante cúbica a trozos de
Hermite de f por x0 , . . . , xn . 4
Observación 3.4.2. En general, la interpolación de Hermite refiere a interpolar una función
f incorporando no solamente sus valores en los nodos (los y0 = f (x0 ), . . . , yn = f (xn ))
sino que también los valores de las derivadas hasta un cierto orden de f en los nodos.
Si queremos que las primeras k derivadas de nuestra interpolante de Hermite a trozos
coincidan con las primeras k derivadas de f , entonces la interpolante debe ser un polinomio
de grado 2k + 1 en cada subintervalo. Por ejemplo, si queremos que nuestra interpolante
tenga iguales derivada primera y derivada segunda en todos los nodos, obtendremos la
llamada interpolante quı́ntica de Hermite, que es un polinomio de grado 5 a trozos. En
este curso nos restringimos a la interpolación cúbica de Hermite, que es la más utilizada
en la práctica. 4
De forma análoga a lo que hicimos en las secciónes 3.3.1 y 3.4.1, podemos estimar el error
de interpolación al usar interpolantes cúbicas a trozos de Hermite. Comenzamos con una
fórmula de representación del error análoga a (3.11).
Teorema 3.4.2 (error de interpolación cúbica de Hermite). Sea f : [a, b] → R una función
de clase C 4 , y consideremos su interpolante cúbica de Hermite p por los puntos a y b, esto
es, p es el único polinomio de grado 3 que cumple p(a) = f (a), p(b) = f (b), p0 (a) = f 0 (a),
p0 (b) = f 0 (b). Entonces, para todo x ∈ [a, b] existe un γx ∈ (a, b) tal que
(x − a)2 (x − b)2 (4)
e(x) := f (x) − p(x) = f (γx ). (3.22)
24
Demostración. Sea x ∈ [a, b]. El resultado es trivial si x = a o x = b, por lo que asumimos
x ∈ (a, b). Consideramos la función auxiliar G : [a, b] → R,
(t − a)2 (t − b)2
G(t) := e(t) − e(x) .
(x − a)2 (x − b)2
Como G(a) = G(b) = G(x) = 0, por el Teorema de Rolle podemos asegurar que G0 tiene al
menos dos raı́ces, una en el intervalo (a, x) y otra en el intervalo (x, b). Pero además, como
86
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Demostración. La demostración es análoga a la del Teorema 3.4.1. Sea x ∈ [xi , xi+1 ] para
algún i ∈ {0, . . . , n − 1}. Para probar (3.23), aplicamos (3.22) en el intervalo [xi , xi+1 ],
tomamos valor absoluto, y como γx ∈ (xi , xi+1 ) podemos acotar
|f (4) (γx )| ≤ kf (4) kL∞ ([xi ,xi+1 ]) .
Luego, para probar (3.24), basta con notar que
(xi+1 − xi )4
(x − xi )2 (xi+1 − x)2 ≤ ∀x ∈ [xi , xi+1 ],
16
y usar que kf (4) kL∞ ([xi ,xi+1 ]) ≤ kf (4) kL∞ ([x0 ,xn ]) y (xi+1 − xi )4 ≤ máxj=0,...,n (xj+1 − xj )4 .
87
Métodos Numéricos Facultad de Ingenierı́a – IMERL
88
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Splines cúbicas
Volvamos a nuestro punto de partida (3.20) y (3.21), pero ahora suponiendo que no
conocemos los valores de d0 , . . . , dn . Esto es, dados n+1 puntos en el plano {(xi , yi )}i=0,...,n
con x0 < x1 < . . . < xn , buscamos una función p cúbica a trozos, derivable en [x0 , xn ], y
tal que p(xi ) = yi para todo i = 0, . . . , n.
La idea al construir splines cúbicas es lograr que la función p sea lo más regular posible.
En el caso de funciones cúbicas a trozos, se toman los {di }i=0,...,n de modo tal que p sea de
clase C 2 , esto es, que tenga derivada segunda continua en el intervalo [x0 , xn ]. Derivando
dos veces (3.21), obtenemos para cada x ∈ [xi , xi+1 ], i = 0, . . . , n − 1,
(6hi − 12s)δi + (6s − 2hi )di+1 + (6s − 4hi )di
p00 (x) = , (3.25)
h2i
donde nuevamente usamos la notación
yi+1 − yi
s = x − xi , hi = xi+1 − xi , δi = .
hi
Como p [xi ,xi+1 ] es un polinomio, para lograr que p sea C 2 basta con imponer que las
derivadas segundas se peguen bien en los nodos interiores x1 , . . . , xn−1 . Usamos (3.25) y
obtenemos, para j = 1, . . . , n − 1,
6δj − 2dj+1 − 4dj
p00 (x+
j ) = (i = j, s = 0 en (3.25)),
hj
−6δj−1 + 4dj + 2dj−1
p00 (x−
j ) = (i = j − 1, s = hi en (3.25)).
hj−1
Por lo tanto, si queremos que la interpolante p sea de clase C 2 , necesitamos que los
coeficientes d0 , . . . , dn satisfagan
6δj − 2dj+1 − 4dj −6δj−1 + 4dj + 2dj−1
= ∀j = 1, . . . , n − 1,
hj hj−1
lo que podemos reescribir como
hj dj−1 + 2(hj−1 + hj )dj + hj−1 dj+1 = 3hj δj−1 + 3hj−1 δj ∀j = 1, . . . , n − 1. (3.26)
Ası́, hemos obtenido un sistema (tridiagonal) de ecuaciones lineales que debemos resolver
para deteminar los valores de d0 , . . . , dn . Recordemos que una vez que tenemos hallados
estos valores, nos basta con utilizar la expresión (3.21) para definir la interpolante spline
cúbica a trozos p.
Sin embargo, notemos que (3.26) pone requerimientos solamente sobre los nodos interiores,
esto es, nos da n − 1 restricciones. En cambio, tenemos n + 1 incógnitas d0 , . . . , dn . Aún
tenemos algo de libertad para imponer dos condiciones adicionales. Aquı́ se abren algunas
posibilidades8 :
8
Referimos a [Sau12, Sección 3.4.2] para más detalles y otras opciones. En ese libro, a la que nosotros
llamamos spline completa se la denomina spline sujetada (clamped).
89
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Prescribir los valores de p0 (x0 ) y p0 (xn ), esto es, tomar d0 y dn como datos, por lo
que ahora (3.26) corresponde a un sistema de n − 1 ecuaciones con n − 1 incógnitas
d1 , . . . , dn−1 . Esta elección se llama spline completa.
Imponer que p00 (x0 ) = 0 y que p00 (xn ) = 0. Esto nos agrega dos ecuaciones adicionales
más a (3.26) y da lugar a un sistema (n + 1) × (n + 1). Hacer esto se conoce como
tomar una spline natural.
Buscar que p tenga un poco más de regularidad cerca de los bordes, esto es, que p
sea de clase C 3 en los intervalos [x0 , x2 ] y [xn−2 , xn ]. Como p es una función cúbica
a trozos, hacer esto es equivalente a usar la misma función cúbica en [x0 , x1 ] que
en [x1 , x2 ] y respectivamente, en [xn−2 , xn−1 ] que en [xn−1 , xn ]. Esto es, estamos
“eliminando” los nodos x1 y xn−1 , y por esta razón a este tipo de interpolante se la
conoce como spline not-a-knot9 .
Ejemplo 3.4.3 (spline not-a-knot). Para fijar ideas, vamos a ejemplificar con la última
elección, la de spline not-a-knot, analizando lo que ocurre en [x0 , x2 ]. Tomando derivada
tercera de p en (3.21) en los intervalos [x0 , x1 ] y [x1 , x2 ], obtenemos dos constantes. Como
queremos que la interpolante tenga derivada tercera continua en x1 , estas dos constantes
deben coincidir. Ası́, obtenemos la identidad
De la misma forma, imponer que p tenga derivada tercera continua en [xn−2 , xn ] da lugar
a la ecuación
−12δn−2 + 6dn−1 + 6dn−2 −12δn−1 + 6dn + 6dn−1
2
= .
hn−2 h2n−1
Por lo tanto, al sistema (3.26) le debemos agregar estas dos ecuaciones y obtenemos un
sistema (n + 1) × (n + 1). 4
90
Métodos Numéricos Facultad de Ingenierı́a – IMERL
En Octave y Matlab hay una función que se llama pchip y que genera una interpolante
“visualmente agradable”10 de los datos. La idea es que la interpolante respete los máximos
y mı́nimos locales, esto es, que no genere oscilaciones más allá de los valores y0 , . . . , yn .
A continuación describimos cómo opera la función pchip para definir los valores de
d1 , . . . , dn y, en consecuencia, construir una interpolante cúbica a trozos. Referimos a
[Mol04, Sección 3.4] para más detalles. Remarcamos, una vez más, que las elecciones que
tomamos aquı́ son arbitrarias y responden a criterios cosméticos.
(2y0 + y1 )δ0 − y0 δ1
d := .
y0 + y1
Luego,
Ejemplo 3.4.4 (spline vs. interpolante “que preserva forma”). Consideremos los datos
de la siguiente tabla:
x 0 1 3 4 5 7
y 1 1 -2 2 3 0
10
Evidentemente, esta es una noción subjetiva.
91
Métodos Numéricos Facultad de Ingenierı́a – IMERL
92
Capı́tulo 4
Ecuaciones no lineales
4.1. Introducción
En este capı́tulo analizamos métodos para aproximar raı́ces de funciones. Para fijar ideas,
comenzamos nuestra discusión con funciones f : R → R, y buscamos x∗ ∈ R tal que
f (x∗ ) = 0. Si f es una función lineal, o más en general si contamos con una forma
de computar exactamente la preimagen de 0, entonces el problema es sencillo. Esto se
corresponde con los métodos directos para resolver sistemas de ecuaciones lineales: por
ejemplo, conocemos una fórmula cerrada para escribir las raı́ces de polinomios cuadráticos
y podrı́amos utilizarla para –al menos en lápiz y papel– determinarlas exactamente. Los
casos en que tenemos una fórmula cerrada para hallar raı́ces son muy particulares, y en
general buscamos una aproximación de x∗ mediante métodos indirectos.
Ası́, tendremos generada una sucesión {xk } que aspiramos que cumpla xk → x∗ con
f (x∗ ) = 0. Como condición de parada, podemos tomar una cierta tolerancia ε > 0 y
computar elementos de esta sucesión hasta que se satisfaga que los errores
ek := xk − x∗ , (4.1)
cumplen |ek | < ε. Este enfoque no parece muy útil en la práctica, ya que para computar
ek debemos conocer x∗ . Una alternativa computable es utilizar los residuos
rk := f (xk )
y parar cuando |rk | < ε. Como ya a esta altura podemos esperar, el hecho de que un
residuo “pequeño” en valor absoluto dé lugar a un error “pequeño” en valor absoluto
depende de un cierto número de condición. Notemos, sin embargo, que tal número de
condición no es igual al de la Definición 1.3.3, ya que nuestro objetivo aquı́ no es evaluar
f (es decir, hallar sus imágenes), sino hallar preimágenes del número 0. Visualmente, la
Figura 4.1 nos muestra que aquı́, si f 0 (x∗ ) ≈ 0, entonces el problema de hallar raı́ces de f
va a estar mal condicionado. Referimos a [Sau12, Sección 1.3.3] para más detalles sobre
el condicionamiento del problema de hallar raı́ces.
93
Métodos Numéricos Facultad de Ingenierı́a – IMERL
2 2
xk xk
0 ∗ 0
x x∗
−2 −2
4
Observación 4.1.1 (velocidad de convergencia). No todos los métodos convergentes tienen
por qué tener un orden o una velocidad de convergencia. En caso de que un método lo
tenga, por lo general conocer el orden nos da mucho mejor idea sobre cómo aproxima a
la solución que conocer su velocidad de convergencia. Cuando p = 1, se suele decir que la
convergencia es lineal. En ese caso, para que sea convergente, es necesario que su velocidad
de convergencia sea β < 1. En caso de que p > 1, no es necesario que sea β < 1 para que
el método sea convergente. 4
94
Métodos Numéricos Facultad de Ingenierı́a – IMERL
En general, una ecuación no lineal puede tener una cantidad arbitraria de soluciones, y el
conjunto de soluciones podrı́a ser arbitrariamente difı́cil de caracterizar matemáticamente.
Asegurar que una ecuación no lineal tiene alguna solución puede ser una tarea difı́cil, pero
el Teorema de Bolzano A.2.1 nos da una situación concreta en la que sı́ podemos lograrlo:
toda función real continua que tenga un cambio de signo en un intervalo [a, b] tiene al
menos una raı́z en dicho intervalo. El Teorema de Bolzano no nos asegura que esa raı́z
sea única –en general, no lo es–, pero sı́ podemos usarlo en forma iterada de modo de
encerrar a alguna raı́z. Esto da lugar al llamado método de bisección.
El método de bisección no es más que una forma sistemática de ensayo y error. Supon-
gamos que I ⊂ R es un intervalo, y que tenemos una función f : I → R continua, con
f (a) < 0, f (b) > 0 para ciertos a, b ∈ I 1 . Podemos considerar el punto medio entre a y b,
m := a+b
2
.
95
Métodos Numéricos Facultad de Ingenierı́a – IMERL
atrapando a x∗ . Vamos a llamar a este tipo de métodos como de encierro2 . Para este
método no tiene sentido definir el error como en (4.1), pero la longitud del intervalo [a, b]
nos da una cota superior para el error. Si usamos como estimación del error la longitud
del intervalo de búsqueda [a, b] y observamos que en cada paso esta longitud se reduce
a la mitad, deducimos que el método de bisección es de primer orden, y tiene
velocidad de convergencia igual a 1/2. En lenguaje más práctico, como log2 (10) ≈ 3,3,
tenemos que cada aproximadamente 3 pasos nuestro intervalo de búsqueda gana un dı́gito
de precisión para aproximar x∗ . 4
Observación 4.2.2 (lento, pero seguro). El método de bisección tiene dos caracterı́sticas
particulares. En primer lugar, tiene la desventaja de que usa información muy cruda de
f : solamente considera su signo y no sus valores. Incluso si f fuese una función lineal
el método no convergirı́a en una sola iteración a menos que la raı́z de f fuese el punto
medio del intervalo inicial de búsqueda. Esta limitación se contrasta con el hecho de que
el método de bisección nos garantiza la convergencia a una raı́z. 4
Observación 4.2.3 (converge, pero, ¿a dónde?). Notemos que si f tiene más de una raı́z
en el intervalo inicial de búsqueda [a, b], entonces puede no ser sencillo caracterizar hacia
qué raı́z converge el método de bisección. 4
96
Métodos Numéricos Facultad de Ingenierı́a – IMERL
3
(a, f (a))
2
x∗
0 mbis mrf
−1 (b, f (b))
0 0,5 1 1,5 2
Figura 4.2: Comparación entre cómo se toma el nuevo punto m para los métodos de
bisección (mbis ) y de la regla falsa (mrf ). A la recta marcada en azul se le llama la secante
al gráfico de f por (a, f (a)), (b, f (b)).
Parece tentador pensar que el método de la regla falsa es de mayor orden que el de
bisección, al estar incorporando información sobre cuánto vale f y no solamente su signo.
Sin embargo, a diferencia de este último, la velocidad de convergencia del método de regla
falsa no siempre es lineal y puede verse afectada por la forma de la función cerca de la
raı́z. En particular, si la función es estrictamente cóncava o convexa en el intervalo, el
método mantendrá una de las aproximaciones fija mientras que la otra se acercará a la
raı́z buscada, lo que puede ralentizar la convergencia. Esta diferencia puede verse en el
siguiente ejemplo.
Ejemplo 4.3.1 (convergencia de regla falsa). Buscamos computar log 3 ≈ 1,0986 a partir
de hallar la raı́z de f (x) = 3e−x − 1. Inicializamos el método de la regla falsa con
a = 0 ⇒ f (a) = 2 > 0,
b = 2 ⇒ f (b) = 3e−2 − 1 ≈ −0,5940 < 0,
ek
n a b m (ver (4.2)) f (m) ek := m − log 3 ek−1
1 0 2 1,5420 −0,3582 0,4434 —
2 0 1,5420 1,3078 −0,1888 0,2092 0,47
3 0 1,3078 1,1950 −0,0919 0,0964 0,46
4 0 1,1950 1,1425 −0,0430 0,0439 0,46
5 0 1,1425 1,1185 −0,0197 0,0199 0,44
6 0 1,1185 1,1076 −0,0089 0,0090 0,45
Hay dos hechos que resultan llamativos. En primer lugar, el punto a permanece fijo en la
iteración, esto es, en todos los pasos tomamos a = 0, pues el punto m computado verifica
f (m) < 0 y por lo tanto al iterado siguiente se toma b ← m. Como a está a una distancia
97
Métodos Numéricos Facultad de Ingenierı́a – IMERL
“grande” de la raı́z log 3, esto implica que si pusiéramos un criterio de parada como en
el Algoritmo 4.1 con tol suficientemente pequeña la iteración no pararı́a nunca4 . En este
sentido, la longitud del intervalo de búsqueda no parece ser un buen indicador de qué tan
cerca está el método de converger. Podrı́a parecer más razonable tomar como indicador
de error la distancia entre m y log 3, como tomamos en la penúltima columna de la tabla;
aquı́ estamos haciendo la “trampa” de que damos por conocido el valor de x∗ que estamos
aproximando. Allı́ surge el segundo hecho llamativo: los errores computados de esta forma
decrecen linealmente, y la velocidad de convergencia parece ser apenas menor que 1/2.
En este ejemplo, no parece claro que el método de la regla falsa se desempeñe mejor que
el de bisección. 4
Ejercicio 4.1. El objetivo de este ejercicio es generalizar el Ejemplo 4.3.1. Consideremos
una función f : R → R de clase C 2 y tal que f 00 > 0.
2. Supongamos que, como en el Ejemplo 4.3.1, existen puntos a < b tales que f (a) > 0,
f (b) < 0 y se inicializa el método de la regla falsa con estos puntos. Demostrar que
en todas las iteraciones se tiene f (m) < 0, y por lo tanto el punto a permanece fijo
a lo largo de la iteración.
3. Análogamente, probar que si existen puntos a < b tales que f (a) < 0, f (b) > 0
y se inicializa el método de la regla falsa con estos puntos, entonces en todas las
iteraciones se tiene f (m) < 0 y por lo tanto b permanece fijo en la iteración.
98
Métodos Numéricos Facultad de Ingenierı́a – IMERL
este método no revisa el signo de f (xk+1 ) para definir a cuál de los puntos xk−1 o xk
descartar. Como este no es un método de encierro, la notación que estamos usando es
un poco distinta; sin embargo, dados puntos a = xk , b = xk−1 , para computar el iterado
siguiente simplemente necesitamos usar (4.2):
99
Métodos Numéricos Facultad de Ingenierı́a – IMERL
L(x) := f (xk ) + sk (x − xk ).
Por otra parte, como L es una interpolante lineal para f en I, podemos utilizar el Teorema
3.3.1 (o la fórmula (3.17) especı́ficamente para una interpolación lineal) para escribir
f 00 (γx∗ ) ∗ f 00 (γk+1 ) k−1 k
L(x∗ ) − f (x∗ ) = − (x − xk−1 )(x∗ − xk ) = − e e , (4.8)
2 2
para algún γk+1 . Combinando (4.7) y (4.8), obtenemos (4.6).
100
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Lema 4.4.3 (Orden de convergencia del método de la secante). Sea f : I → R una función
de clase C 2 , y x∗ un punto interior a I y tal que f (x∗ ) = 0, f 0 (x∗ ) 6= 0.
√
Entonces, cuando
1+ 5
el método de la secante converge con orden, éste es al menos Φ := 2 .
Ahora, como sabemos que el método converge con orden, podemos asumir que para k
suficientemente grande se cumplen |ek+1 | ≈ C|ek |p y |ek | ≈ C|ek−1 |p , donde p es el orden
de convergencia y C es la velocidad de convergenica del método. Reemplazando arriba,
tenemos
1 p+1
C|ek |p . M C|ek−1 |p+1 ⇒ |ek | . M p |ek−1 | p .
1 p+1
Como tenemos |ek | ≈ C|ek−1 |p y |ek | . M p |ek−1 | p , deducimos que el orden p debe
cumplir
p+1
p= , o equivalentemente p2 − p − 1 = 0.
p
√
1+ 5
Concluimos que el orden de convergencia es p = Φ = 2
.
Observación 4.4.2 (cuenca de convergencia). El Teorema 4.4.4 nos dice que, bajo condi-
ciones adecuadas, si el método de la secante es inicializado lo suficientemente cerca de
la raı́z entonces es convergente. Notemos que a partir de la demostración se deduce una
forma para estimar un valor de δ > 0 que garantice la convergencia: basta con tomar δ
de modo que la derivada f 0 varı́e lo suficientemente poco (en el sentido de (4.5)) en un
entorno de centro x∗ y radio δ. De todos modos, estimar el valor de δ es un tema delicado
en general, ya que si no conocemos x∗ , parece difı́cil que podamos conocer f 0 (x∗ ) con
exactitud. 4
Observación 4.4.3 (velocidad de convergencia). La Proposición 4.4.2, junto con la demos-
tración del Lema 4.4.3, permiten incluso estimar la velocidad de convergencia del método.
101
Métodos Numéricos Facultad de Ingenierı́a – IMERL
En efecto, si ahora ya asumimos que el método converge con orden Φ y volvemos a (4.6),
podemos afinar nuestra cota para M y aproximar
f 00 (x∗ )
f 00 (γk+1 ) ≈ f 00 (x∗ ), sk ≈ f 0 (x∗ ) ⇒M ≈ .
2f 0 (x∗ )
Como la velocidad de convergencia que estimamos en la demostración del Lema 4.4.3 era
1
C = M Φ , podemos estimar la velocidad del método como
1
f 00 (x∗ ) Φ
C≈ .
2f 0 (x∗ )
4
Ilustramos el Teorema 4.4.4 con el siguiente ejemplo.
Ejemplo 4.4.1 (convergencia de secante). Repetimos la consigna del Ejemplo 4.3.1, pero
ahora usamos el método de la secante inicializado en los mismos puntos x0 = 0, x1 = 2.
Computamos algunos iterados en la siguiente tabla.
|ek+1 |
k xk−1 xk xk+1 (ver (4.4)) ek+1 := xk+1 − log 3 |ek |Φ
1 0 2 1,5420 0,4434 0,5245
2 2 1,5420 0,8465 −0,2521 0,9398
3 1,5420 0,8465 1,1557 0,0571 0,5311
4 0,8465 1,1557 1,1056 0,0070 0,7144
5 1,1557 1,1056 1,0984 −0,0002 0,6221
6 1,1056 1,0984 1,0986 7 × 10−7 0,6712
La última columna nos indica que el orden de convergencia observado es Φ, con una veloci-
dad de convergencia alrededor de 0,67. Por la Observación 4.4.3, como en este ejemplo te-
nemos f 0 (x∗ ) = −1, f 00 (x∗ ) = 1, esperamos convergencia con velocidad 2−1/Φ ≈ 0,65. 4
102
Métodos Numéricos Facultad de Ingenierı́a – IMERL
2
(xk , f (xk ))
1
xk+1
0
−1
0 0,5 1 1,5 2
Como este polinomio es lineal, resulta trivial hallar su raı́z y ası́ queda definido xk+1 de
la misma forma que en (4.9). 4
Observación 4.5.2 (comparación con secante). El método de Newton-Raphson no solo
utiliza los valores de f , sino que también incorpora información sobre f 0 . Por lo tanto,
es de esperar que presente una convergencia de mayor orden que el de la secante. Sin
embargo, requiere que tengamos acceso a f 0 , lo que no siempre tiene por qué ocurrir en
la práctica. 4
103
Métodos Numéricos Facultad de Ingenierı́a – IMERL
∗ k 0 k ∗ f 00 (θk ) ∗
k
0 = f (x ) = f (x ) + f (x )(x − x ) + (x − xk )2 ,
2
para algún θk entre x∗ y xk . Como f 0 es continua (pues f ∈ C 2 ), f 0 (x∗ ) 6= 0, y el método
es convergente, existe k0 tal que f 0 (xk ) 6= 0 para todo k ≥ k0 . Combinamos la definición
(4.9) con el desarrollo de Taylor y obtenemos
104
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(xk )γ
k+1 k k 1
x =x − =x 1− .
γ(xk )γ−1 γ
ek+1 1
k
=1− ,
e γ
ek+1 1
k
=1− .
e γ
1
Como 0 < γ < 1, tenemos 1 − γ
< 0. De hecho, si 1/2 < γ < 1 aún tenemos
1
1− γ
< 1 y podemos garantizar que el método converge linealmente. Si 0 < γ <
1
1/2, entonces 1 − γ
> 1 y tendremos |ek | → ∞. Finalmente, en el caso lı́mite
γ = 1/2, la iteración de Newton-Raphson permanecerá “rebotando” entre x0 y −x0 .
105
Métodos Numéricos Facultad de Ingenierı́a – IMERL
¡Una perturbación del orden de 10−3 en el iterado inicial afecta a qué raı́z converge el
método! 4
106
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Nos referiremos a un método que se pueda escribir de esta forma como un método
iterativo general.
Naturalmente, dada una función f : R → R, existen infinitas formas posibles de elegir
g de modo que se cumpla (4.10). Una de estas formas es g(x) := x − ff0(x)(x)
, que, como
observamos, corresponde al método de Newton-Raphson y muestra que éste es un método
iterativo general. Otra posibilidad para construir un MIG a partir de una función f , es
tomar, por ejemplo, g(x) := x − f (x).
La Observación 4.5.3 muestra que no siempre podemos esperar que la iteración (MIG)
sea convergente. Sin embargo, el Teorema 4.5.1 sugiere que la derivada de la función
f (y por lo tanto de la función g para nuestro MIG) debe cumplir un rol fundamental
en la convergencia de estos métodos. Como muestra el Corolario 4.7.2 (ver abajo), para
que (MIG) sea convergente basta con asegurar que |g 0 | sea menor que 1 a lo largo de la
iteración.
x2 − x − 2 = 0 ⇔ x = x2 − 2 g1 (x) := x2 − 2,
√
x2 − x − 2 = 0 ⇔ x2 = x + 2 g2 (x) := x + 2,
2 2
x2 − x − 2 = 0 ⇔ x(x − 1) = 2 ⇔ x = 1 + g3 (x) := 1 + ,
x x
2
x +2
x2 − x − 2 = 0 ⇔ x2 − x = 2 ⇔ x(2x − 1) = x2 + 2 g4 (x) := .
2x − 1
107
Métodos Numéricos Facultad de Ingenierı́a – IMERL
4.7.1. Convergencia
La siguiente definición nos permite dar una condición compacta para asegurar la conver-
gencia de un método iterativo general.
b) Probar que si g es una función derivable en I, y existe un λ < 1 tal que |g 0 (x)| < λ
para todo x ∈ I, entonces g es contractiva en I.
El siguiente es un resultado muy general, que en nuestro contexto nos permite asegurar
la convergencia de (MIG).
Teorema 4.7.1 (de punto fijo). Sea I ⊂ R un intervalo cerrado, y sea g : I → I una
función contractiva. Entonces, existe un único x∗ ∈ I tal que g(x∗ ) = x∗ . Más aún, para
todo x0 ∈ I, si se define la sucesión {xk }k∈N mediante (MIG), se tiene
lı́m xk = x∗ .
k→∞
108
Métodos Numéricos Facultad de Ingenierı́a – IMERL
109
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ya tenemos los ingredientes necesarios para dar condiciones que nos permitan asegurar la
convergencia de (MIG).
xk ∈ Iδ ∀k, lı́m xk = x∗ .
k→∞
Demostración. Sea x0 ∈ Iδ . Nos basta con hallar un intervalo cerrado I˜δ ⊂ Iδ que contenga
a x0 y a x∗ . En efecto, supongamos I˜δ = [x∗ − δ̃, x∗ + δ̃] para algún 0 < δ̃ < δ. Como g es
contractiva en Iδ , es fácil verificar que también lo es en I˜δ . Por lo tanto, para todo y ∈ I˜δ ,
tenemos
y concluimos g(y) ∈ I˜δ para todo y ∈ I˜δ . Como g|I˜δ verifica las hipótesis del Teorema
4.7.1, concluimos que (MIG) converge al único punto fijo de g en I˜δ , que es x∗ .
A continuación, damos un caso particular del Corolario 4.7.2 que es de utilidad práctica.
Demostración. Al ser g 0 una función continua y tener |g 0 (x∗ )| < 1, existen λ < 1 y δ > 0
tal que si |x − x∗ | < δ, entonces |g 0 (x)| < λ. Por el Ejercicio 4.2, parte b), deducimos que
g es contractiva en (x∗ − δ, x∗ + δ) y podemos aplicar el Corolario 4.7.2.
Ejemplo 4.7.2 (distintos MIGs, revisitado). Volvemos al Ejemplo 4.7.1, en el que bus-
camos la raı́z x∗ = 2. Tenemos
g10 (x) = 2x − 1 ⇒ g10 (x∗ ) = 3,
g20 (x) = 2√x+2
1
⇒ g20 (x∗ ) = 14 ,
g30 (x) = − x22 ⇒ g30 (x∗ ) = − 12 ,
2 +2)
g40 (x) = 2x(2x−1)−2(x
(2x−1)2
⇒ g40 (x∗ ) = 0.
Esto muestra que la iteración usando g10 no es convergente, mientras las otras tres sı́ lo son,
siempre y cuando las inicialicemos lo suficientemente cerca de x∗ . Además, observamos
que cuanto menor sea |gi0 (x∗ )| más rápida parece ser la convergencia. 4
110
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Demostración. Ya sabemos que (4.9) es de la forma (MIG), tomando g(x) = x − ff0(x) (x)
.
0 ∗
Nos basta con verificar que se verifica |g (x )| < 1 y aplicar el Corolario 4.7.3. Derivamos
g, y obtenemos
0 f 0 (x)f 0 (x) − f 00 (x)f (x) f 00 (x)f (x)
g (x) = 1 − = .
f 0 (x)2 f 0 (x)2
Por lo tanto, g 0 (x∗ ) = 0.
para algún θx entre x∗ y x. Evaluamos este desarrollo en xk , notando que xk+1 = g(xk ),
para obtener
g (p) (θk ) k
xk+1 = g(xk ) = x∗ + (x − x∗ )p .
p!
Reordenando, esta relación se puede escribir en función de los errores como
|g(p) (x∗ )|
lo que demuestra que la convergencia es con orden p y velocidad p!
.
111
Métodos Numéricos Facultad de Ingenierı́a – IMERL
|g(p) (x∗ )|
(⇒) Supongamos ahora que (MIG) converge con orden p y velocidad p!
, lo que
significa que
ek+1 g (p) (x∗ )
lı́m p = 6= 0.
k→∞ |ek | p!
Esto implica que g (p) (x∗ ) 6= 0. Veamos ahora que g (i) (x∗ ) = 0, ∀ i = 1, . . . , p − 1. Supon-
gamos por absurdo que alguna de estas derivadas no es nula, y sea j el orden de la primer
derivada no nula en x∗ , esto es,
112
Métodos Numéricos Facultad de Ingenierı́a – IMERL
¡Nuestro control sobre el error es mucho más pobre! Naturalmente, este razona-
miento se puede generalizar: si f tiene sus primeras p − 1 derivadas nulas en
x∗ y f (p) (x∗ ) √
6= 0, entonces si |f (xk )| < ε solamente podremos afirmar que
|xk − x∗ | = O( p ε).
Usar el número de iteraciones. Esta puede ser una importante red de seguridad
en caso de que el método no esté convergiendo o lo haga muy lentamente. Podemos
fijar que el algoritmo se detenga en caso de que el número de iteraciones llegue a
una cantidad máxima fijada de antemano. De este modo, evitamos que el algoritmo
entre en un ciclo sin fin, pero naturalmente no tenemos ningún tipo de información
relevante sobre x∗ en caso de que el algoritmo se detenga usando este criterio.
113
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Consideramos f : Rn → Rn tal que f (x) = [f1 (x), f2 (x), . . . , fn (x)]t . En forma compacta,
resolver (4.12) puede expresarse equivalentemente como
Para simplificar la discusión, de aquı́ en más asumimos que g es tan regular como sea
necesario. No vamos a profundizar en aspectos técnicos de este tipo de métodos, pero
con los conocimientos que hemos adquirido en la Sección 4.7, podemos generar algunas
intuiciones. Por una parte, en el caso de métodos iterativos generales para ecuaciones no
lineales en R tenemos que la derivada de la función de iteración juega un rol prepon-
derante para establecer la convergencia, asegurando la misma en caso de que la función
sea contractiva. Para funciones de varias variables, la matriz jacobiana cumple el rol de
la derivada, y por lo tanto es de esperar que deba asignarse una condición que asegure
contractividad para la matriz jacobiana.
Por otro lado, de nuestros conocimientos de Geometrı́a y Álgebra Lineal 2 sabemos que
la intuición de “contraer direcciones” tiene que ver con la noción de vectores y valores
propios. La siguiente noción nos permite tratar esta idea en forma compacta.
114
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Definición 4.8.1 (radio espectral). Sea A ∈ Mn (R) una matriz cuyos valores propios
(complejos) escribimos como λ1 , . . . , λn . El radio espectral de A, que denotamos por
ρ(A), se define como
ρ(A) = máx |λi |.
1≤i≤n
Intuitivamente, la Definición 4.8.1 nos dice que el hecho de que una matriz A ∈ Mn (R)
“contraiga” todas las direcciones es equivalente a que ρ(A) ≤ 16 .
Combinando la definición de radio espectral con la matriz jacobiana, y en analogı́a con
el Teorema 4.7.1, es de esperar que podamos garantizar que (4.14) converja a x∗
siempre y cuando la matriz jacobiana de g tenga radio espectral menor que
1 en un entorno de x∗ , y de hecho, cuanto menor sea ρ(Jg (x∗ )) más rápida será la
convergencia.
cuya solución x∗ = [0, 1]t buscamos aproximar, y analizamos dos alternativas usando
(4.14). Consideramos g1 , g2 : R2 → R2 , tales que
" #
−2x + 2
2x + 2x − 2 √ 22
g1 (x1 , x2 ) = 2 1 2 2 , g2 (x1 , x2 ) = 4−x1 .
x1 + 4x2 + x2 − 4
2
Notemos que ambas alternativas satisfacen (4.13), pero sus jacobianos son
2 2 ∗ 2 2
Jg1 (x1 , x2 ) = ⇒ Jg1 (x ) = ,
2x1 8x2 + 1 0 9
" #
0 −2
∗ 0 −2
Jg2 (x1 , x2 ) = √−x1 0 ⇒ Jg2 (x ) = 0 0 .
2 4−x2 1
Es sencillo verificar que los valores propios de Jg1 (x∗ ) son 9 y 2, mientras que Jg2 (x∗ )
solamente tiene valor propio 0 (que es doble). Deducimos que ρ(Jg1 (x∗ )) = 9 y que
ρ(Jg2 (x∗ )) = 0, por lo que esperamos que la iteración usando g1 no converja y que la
iteración usando g2 sı́ lo haga.
Implementamos en Octave estas iteraciones, comenzando con x0 = [0,5, 0,5]t . La siguiente
tabla muestra los primeros iterados e indica que nuestro razonamiento parece ser el co-
rrecto. Para no cargar la notación, usamos un punto y coma ; para marcar un cambio de
fila.
6
Es claro por definición que ρ(A) ≥ 0 para toda matriz A ∈ Mn (R).
115
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Observación 4.8.1 (criterios de parada). Los criterios de parada que imponemos sobre los
métodos iterativos en Rn son los mismos que comentamos en la Sección 4.7.3, con la única
diferencia de que donde dice “valor absoluto” allı́, aquı́ debemos usar alguna norma en
Rn7 . 4
De este modo, dado xk , para para calcular el siguiente iterado xk+1 podemos
116
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Observación 4.8.2 (invertibilidad de Jf ). Para que los iterados del método de Newton-
Raphson estén bien definidos, es necesario que las matrices jacobianas Jf (xk ) sean inver-
tibles. Esto es análogo a la versión (4.9) en R, en la que necesitamos que f 0 (xk ) 6= 0. 4
Ejemplo 4.8.2 (Newton-Raphson en R2 ). Volvemos al problema del Ejemplo 4.8.1, al
que ahora abordamos usando el método de Newton-Raphson en R2 . Para computar los
iterados necesitamos la matriz jacobiana de f , que es
1 2
Jf (x1 , x2 ) = .
2x1 8x2
y por lo tanto
1 0 −5/4
1
x =x +d = .
13/8
Siguiendo de este modo, implementando el Algoritmo 4.2 en Octave obtenemos los iterados
1 −1,2500 2 −0,3472 3 −0,0447 4 −0,0010
x = , x = , x = , x = .
1,6250 1,1736 1,0224 1,0005
117
Métodos Numéricos Facultad de Ingenierı́a – IMERL
118
Métodos Numéricos Facultad de Ingenierı́a – IMERL
secante son métodos de cuasi-Newton que logran obtener órdenes de convergencia mayor
a 1, y toman B k de modo que se cumpla
119
Métodos Numéricos Facultad de Ingenierı́a – IMERL
120
Capı́tulo 5
Mı́nimos cuadrados
yi ≈ f (ti ) := x1 ti + x2 , i = 1, . . . , m.
x1 t2i + x2 ti yi + x3 yi2 + x4 ti + x5 yi + x6 ≈ 0, i = 1, . . . , m,
121
Métodos Numéricos Facultad de Ingenierı́a – IMERL
y y
t t
122
Métodos Numéricos Facultad de Ingenierı́a – IMERL
y tomamos
t−µ1 2
t−µr 2
+ . . . + x e−( σr ) .
−
f (t) = x e σ1
1 r
En las primeras secciones del capı́tulo, nos vamos a concentrar en problemas lineales:
dadas n funciones linealmente independientes φ1 (t), . . . , φn (t), buscamos una función f
de la forma
Xn
f (t) = xj φj (t).
j=1
Notemos que la matriz de diseño almacena las evaluaciones de todas las funciones de base
φ1 , . . . , φn en todos los puntos t1 , . . . , tm . Por lo tanto, si tenemos un problema lineal y
queremos f (ti ) ≈ yi para todo i = 1, . . . , m, estamos buscando un vector de parámetros
x = [x1 , . . . , xn ]t tal que
x1 φ1 (t1 ) + x2 φ2 (t1 ) + . . . + xn φn (t1 ) = y1
x1 φ1 (t2 ) + x2 φ2 (t2 ) + . . . + xn φn (t2 ) = y2
.. .. . (5.1)
. .
x φ (t ) + x φ (t ) + . . . + x φ (t ) = y
1 1 m 2 2 m n n m m
En forma compacta, este sistema se puede escribir como Ax = y, donde A ∈ Mm×n (R) es
la matriz de diseño, e y = [y1 , . . . , ym ]t . A diferencia de lo que estudiamos en los capı́tulos
2 y 3, este no es un sistema cuadrado; como m > n, en general esperamos que este sistema
sea incompatible2 .
2
Notemos que, si m = n, el problema (5.1) es un problema de interpolación, aunque no necesariamente
de interpolación polinomial.
123
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Dado que en general no tenemos esperanza de resolver el sistema (5.1) en forma exac-
ta, necesitamos dar una noción de solución con la que podamos determinar valores
“satisfactorios” de los parámetros. Dado un vector x ∈ Rn , consideremos el residuo
r := Ax − y ∈ Rm . Parece natural buscar un vector de parámetros x ∈ Rn que minimice
alguna norma del residuo. Para esto, puede ser útil recordar las normas `p que intro-
dujimos en la Sección 2.5.1: apuntar a minimizar las normas `1 o `∞ del residuo tiene
aplicaciones en la práctica y pueden entenderse como problemas de programación lineal.
En este capı́tulo nos vamos a concentrar en el problema de minimizar la norma euclı́dea
(`2 ) del residuo.
Definición 5.1.2 (mı́nimos cuadrados). Dados A ∈ Mm×n (R), y ∈ Rm , decimos que
un vector x∗ ∈ Rn es solución del sistema Ax = y en el sentido de mı́nimos
cuadrados si
x∗ = arg mı́nn kAx − yk22 .
x∈R
4
Observación 5.1.1 (cuadrado de la norma). Como seguramente el lector ya se haya encon-
trado en algún curso de cálculo, es fácil comprobar que el punto en el que se minimiza la
norma euclı́dea es el mismo que en el que se minimiza el cuadrado de la norma euclı́dea.
En la definición anterior usamos el cuadrado de la norma euclı́dea simplemente porque es
más cómodo para trabajar. 4
Observación 5.1.2 (mı́nimos cuadrados con pesos). A veces puede ocurrir que queramos
darle más importancia al ajuste de ciertos puntos que a otros. Por ejemplo, esto puede
ocurrir si los {(ti , yi )} corresponden a observaciones experimentales y algunas de estas
observaciones fueron tomadas en condiciones más favorables, o con mejor equipamiento.
En ese caso, se puede tomar w1 , . . . , wm > 0 y apuntar a minimizar la función
m
X
Φw (x) := wi ri2 .
i=1
124
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Nos resta verificar que tal punto crı́tico efectivamente es un mı́nimo de Φ. Sea x∗ tal que
At Ax∗ = At y y sea z ∈ Rn arbitrario. Usamos la condición de punto crı́tico para escribir
Φ(z) − Φ(x∗ ) = zt At Az − 2zt At y + yt y − x∗t At Ax∗ − 2x∗t At y + yt y
125
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Observación 5.2.1 (unicidad). Remarcamos que el Teorema 5.2.1 no indica que el problema
de mı́nimos cuadrados tenga solución única. De hecho, esto solamente ocurre si A es de
rango completo (esto es, rg(A) = n): en ese caso, At A ∈ Mn (R) también tendrı́a que
tener rango igual a n, lo que implica que debe ser invertible y por lo tanto las ecuaciones
normales tienen solución única. En este caso tenemos que At A es definida positiva e incluso
podrı́amos usar el criterio de la Hessiana para clasificar el único punto crı́tico de Φ: basta
con notar que D2 Φ(x) = At A para todo x ∈ Rn .
En caso de que A sea de rango deficiente (esto es, rg(A) < n), entonces At A no es inverti-
ble, y las ecuaciones normales constituyen un sistema indeterminado. En la Sección 6.2.2
del Capı́tulo 6 analizamos este tipo de problemas. Por el momento, solamente comenta-
mos que en este caso tenemos que At A es semidefinida positiva y por lo tanto el uso del
criterio de la Hessiana en la demostración del Teorema 5.2.1 requerirı́a algo de cuidado
adicional. En cambio, enfatizamos que el argumento que hicimos en (5.3) para probar que
Φ alcanza un mı́nimo en x∗ es válido aún en este caso.
4
Observación 5.2.2 (ortogonalidad). La condición At Ax∗ = At y es equivalente a que Ax∗ −
y sea ortogonal al espacio de columnas de A. En efecto, las ecuaciones normales se pueden
escribir como
At (Ax∗ − y) = 0.
Esto implica que, para toda fila de At (esto es, para toda columna A(i) de A), se tiene
A(i) · (Ax∗ − y) = 0. 4
Observación 5.2.3 (mı́nimos cuadrados como proyección ortogonal). Supongamos que m >
n y que rg(A) = n. El espacio de columnas de A (col(A)) es un subespacio de Rm
de dimensión n. Descomponemos Rm como la suma directa entre este subespacio y su
complemento ortogonal,
Rm = col(A) ⊕ col(A)⊥ .
La observación anterior nos asegura que y − Ax∗ ∈ col(A)⊥ , mientras que es evidente que
Ax∗ ∈ col(A). Por lo tanto, la solución de mı́nimos cuadrados es el único vector x∗ ∈ Rn
que nos permite descomponer
Ax∗ + y − Ax∗ .
y = |{z}
| {z }
∈col(A) ∈col(A)⊥
Esto implica que Ax∗ es la proyección ortogonal (ver Definición A.4.2) del
vector y sobre col(A). 4
126
Métodos Numéricos Facultad de Ingenierı́a – IMERL
t y
0 -1
1 1
2 2
3 9/2
Queremos hacer un ajuste lineal a los mismos, esto es, escribir y(t) = x1 t + x2 con x1 , x2
a determinar. Usamos las funciones de base φ1 (t) = t, φ2 (t) = 1, y escribimos nuestro
problema como Ax = y, donde
0 1 −1
1 1 1 x
A= , y = , x = 1 .
(5.4)
2 1 2 x2
3 1 9/2
La solución de este último sistema es x∗ = [7/4, −1]t . La Figura 5.2 muestra los datos y
el gráfico de la función lineal y(t) = 7t/4 − 1 que obtuvimos. 4
Observación 5.2.4 (condicionamiento). Las ecuaciones normales son una herramienta po-
derosa para trabajar en la teorı́a de problemas de mı́nimos cuadrados. A primera vista,
incluso parecen interesantes en la práctica: estamos convirtiendo el problema de resolver
el sistema Ax = y, con A ∈ Mm×n (R) en el de resolver un sistema cuadrado cuya matriz
es At A ∈ Mn (R). En caso de que m n, esto puede ser un buen beneficio respecto
al tamaño del sistema resultante. Sin embargo, las ecuaciones normales no siempre
127
Métodos Numéricos Facultad de Ingenierı́a – IMERL
son una buena opción práctica, y esto se debe a que pueden presentar pro-
blemas de condicionamiento. Para sistemas pequeños, esto puede no ser demasiado
problemático (a menos que tengamos un ejemplo patológico como el que mostramos a
continuación).
Sin embargo, dado que los datos de entrada en general están sujetos a errores (en particu-
lar, si corresponden a observaciones empı́ricas o datos experimentales), en problemas de
mı́nimos cuadrados es especialmente importante reducir la magnificación del error. Por
ello, en la Sección 5.3 y el Capı́tulo 6 vamos a estudiar métodos para resolver problemas
de mı́nimos cuadrados que no requieran el uso de la matriz At A. 4
Ejemplo 5.2.2 (condicionamiento de ecuaciones normales). Dado δ > 0, consideremos
1 1
A = δ 0 ∈ M3×2 (R).
0 δ
5.3. Descomposición QR
Aquı́ analizamos una alternativa para resolver problemas de mı́nimos cuadrados que no
implique utilizar las ecuaciones normales. Consideramos la ecuación Ax = y con A ∈
Mm×n (R), m > n. Como punto de partida, retomamos la Observación 5.2.3: hallar la
solución del sistema en el sentido de los mı́nimos cuadrados es equivalente a computar la
proyección ortogonal de y sobre el espacio de columnas de A. Esto motiva a los llamados
métodos de ortogonalización para resolver nuestro problema.
El lector debe haberse encontrado con el concepto de matriz ortogonal en los cursos de
Geometrı́a y Álgebra Lineal. Aquı́ hacemos un breve repaso de propiedades clave para
nuestro análisis.
128
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Definición 5.3.1 (matriz ortogonal). Sea Q ∈ Mn (R). Decimos que Q es una matriz
ortogonal si sus columnas {Q(1) , . . . , Q(n) } forman una base ortonormal de Rn (recordar
la Definición A.4.1). 4
Observación 5.3.1 (propiedades de matrices ortogonales). Mencionamos algunas propie-
dades que vamos a utilizar sobre matrices ortogonales.
1. Las transformaciones lineales dadas por matrices ortogonales preservan norma. Esto
es, si Q ∈ Mn (R) es ortogonal, entonces se tiene
kQxk2 = kxk2 ∀x ∈ Rn .
4
Teorema 5.3.1 (descomposición QR). Dada A ∈ Mm×n (R) con m ≥ n tal que rg(A) =
n, existen
tales que
A = QR. (5.5)
129
Métodos Numéricos Facultad de Ingenierı́a – IMERL
130
Métodos Numéricos Facultad de Ingenierı́a – IMERL
donde z(1:n) ∈ Rn y z(n+1:m) ∈ Rm−n , esto es, z(1:n) tiene las primeras n coordenadas de z
y z(n+1:m) las últimas m − n. Es simple verificar que4
kzk22 = kz(1:n) k22 + kz(n+1:m) k22 .
En particular, observemos que al aplicar esta forma de descomponer vectores nos da
(Rx)(1:n) = R(1:n) x y (Rx)(n+1:m) = O(m−n)×n x = 0 y por lo tanto, volviendo a (5.6),
tenemos
krk22 = kRx − Qt yk22 = kR(1:n) x − (Qt y)(1:n) k22 + k(Qt y)(n+1:m) k22 . (5.7)
Esto quiere decir que el cuadrado de la norma euclı́dea del residuo –cantidad que queremos
minimizar– se puede escribir como la suma de dos términos: el primero depende del vector
de parámetros x, pero el segundo no. Por lo tanto, la solución del problema de mı́nimos
cuadrados es aquel x∗ ∈ Rn que logre que k(Rx)(1:n) − (Qt y)(1:n) k22 sea lo menor posible.
Sin embargo, el sistema
R(1:n) x = (Qt y)(1:n)
es un sistema cuadrado y triangular superior: lo podemos resolver usando una sustitución
hacia atrás (Sección 2.2.1) sin pasar por ningún problema de condicionamiento. Resolver
este sistema nos da la solución x∗ ∈ Rn .
En resumen, dada la factorización A = QR de una matriz de rango completo, para hallar
x∗ = arg mı́nn kAx − yk22 ,
x∈R
computamos el vector Qt y;
resolvemos el sistema Rx = (Qt y)(1:n) usando sustitución hacia atrás para hallar la
solución de mı́nimos cuadrados x∗ ;
si nos interesa el valor mı́nimo de la norma euclı́dea del residuo, ésta es igual a
k(Qt y)(n+1:m) k2 .
Ejemplo 5.3.1 (ajuste lineal con QR). Volvemos a hacer el Ejemplo 5.2.1, pero esta
vez aplicando el método que recién describimos. Usamos el comando [Q,R] = qr(A) en
Octave5 para obtener
0 −0,8367 −0,3858 −0,3888 −3,7417 −1,6036
−0,2673 −0,4781 0,2246 0,8060 0 −1,1952
Q= −0,5345 −0,1195 0,7082 −0,4455
, R = .
0 0
−0,8018 0,2390 −0,5470 0,0283 0 0
Por lo tanto, podemos hallar la solución de nuestro problema con los comandos
4
Con la sutileza de que las tres normas que aparecen en la igualdad son en espacios distintos: la
primera es en Rm , la segunda en Rn , y la tercera en Rm−n .
5
En la Sección 5.3.2 estudiamos cómo computar esta factorización en forma estable.
131
Métodos Numéricos Facultad de Ingenierı́a – IMERL
>> z = Q’*y;
>> x = R(1:2, 1:2) \ z(1:2);
Observación 5.3.3 (problemas con rango deficiente). Si rg(A) < n en el Teorema 5.3.1,
aún ası́ existe la factorización QR, pero naturalmente en ese caso no podemos esperar que
R sea invertible. De acuerdo a lo que comentamos en la Observación 5.2.1, un problema
de mı́nimos cuadrados que involucre a la matriz A no tiene solución única. En la práctica,
el rango de una matriz no es algo tan claro; una matriz que con aritmética real es de
rango completo, no tiene por qué serlo en aritmética de punto flotante. Si un problema de
mı́nimos cuadrados está cerca de ser de rango deficiente, en la práctica se pueden utilizar
variantes al método de computar la factorización QR como realizar pivoteo por columnas.
No vamos a tratar ese tipo de métodos en este curso, y para problemas de rango deficiente
o casi de rango deficiente proponemos aplicar la factorización SVD que describiremos en
el Capı́tulo 6. 4
H := I − ρuut , (5.8)
donde ρ := 2/kuk22 . 4
En forma un tanto vaga, vamos a usar la expresión reflexión de Householder para referirnos
tanto a la transformación lineal T : Rm → Rm dada por T x = Hx como a su matriz
132
Métodos Numéricos Facultad de Ingenierı́a – IMERL
asociada H. Las reflexiones de Householder son, entonces, matrices que se pueden escribir
uut
como la resta de la matriz identidad en Mm (R) menos la matriz 2 kuk 2 , que es de rango
2
6 u u
uno , y que también se suele escribir como 2 kuk2 ⊗ kuk2 .
Observación 5.3.4. En este punto, puede ser instructivo computar las imágenes de algunos
vectores por una transformación de Householder. Sea H la transformación de Householder
asociada a un cierto vector u ∈ Rm . Podemos calcular imágenes usando la definición (5.8).
Por un lado, vectores de la forma v = αu con α ∈ R se transforman en
Hv = (I − ρuut )v = v − ρuut v = v.
133
Métodos Numéricos Facultad de Ingenierı́a – IMERL
kHxk2 = kxk2 ∀x ∈ Rm .
−PU (x)
u x − PU (x)
Hx = x − 2PU (x)
Figura 5.3: La imagen de un vector x mediante una reflexión de Householder con vector
u es la simetrı́a de x respecto al subespacio [u]⊥ .
134
Métodos Numéricos Facultad de Ingenierı́a – IMERL
135
Métodos Numéricos Facultad de Ingenierı́a – IMERL
De acuerdo a nuestra discusión anterior, basta con que H1 sea la transformación de Hou-
seholder asociada al vector u := A(1) ± kA(1) k2 e1 , donde el signo del segundo sumando
se suele tomar igual que el del primer elemento de A(1) . Luego, al hacer esta primera
operación, tendremos
±kA(1) k2 ∗ ... ∗ ±kA(1) k2 ∗ . . . ∗
0 ∗ . . . ∗ 0
A2 := H1 A = = ∈ Mm×n (R).
.. .. Ã2
. ∗ . . . ∗ .
0 ∗ ... ∗ 0
Los ∗ arriba denotan elementos que no necesariamente son nulos. Observemos que la
primera columna de A2 solamente tiene un elemento no nulo en su primer entrada, por lo
que es la primera columna de una matriz triangular superior. Sea Ã2 ∈ M(m−1)×(n−1) (R)
la matriz que se obtiene al quitar la primera fila y la primera columna de A2 .
A continuación, consideramos una transformación de Householder H̃2 ∈ Mm−1 (R) que
mande la primera columna de Ã2 en el primer vector de la base canónica de Rm−1 ,
(1) (1)
H̃2 Ã2 = ±kÃ2 k2 e1 ,
con Ã3 ∈ M(m−2)×(n−2) (R). Notemos que las dos primeras columnas en A3 se correspon-
den a las primeras dos columnas de una matriz triangular superior. Siguiendo con este
procedimiento, luego de n pasos se obtiene (5.9). Concretamente, en el paso k-ésimo,
136
Métodos Numéricos Facultad de Ingenierı́a – IMERL
tenemos
±kA(1) k2 . . . ∗ ∗ ... ∗
...
0 ∗ ∗ ... ∗
(1)
0 . . . kÃk−1 k2 ∗ . . . ∗
Ak := ,
0 ... 0
.. ..
Ãk
. .
0 ... 0
y tomamos la matriz ortogonal y simétrica
Ik−1 O(k−1)×(m−k+1)
Hk := ,
O(k−1)×(m−k+1) H̃k
Ax = y ⇒ Hn . . . H1 Ax = Hn . . . H1 y =: z ⇒ Rx = z.
La matriz R ∈ Mm×n (R) resultante es triangular superior, por lo que estamos nuevamente
en una situación como en la Sección 5.3.1: para hallar x, nos basta con restringirnos a
las primeras n filas del sistema anterior y hacer una sustitución hacia atrás. La norma
euclı́dea del residuo mı́nimo es igual a la norma euclı́dea del vector formado por las últimas
n − m + 1 entradas de z.
137
Métodos Numéricos Facultad de Ingenierı́a – IMERL
√
u := A(1) + 14e1 , esto es,
0 −0,2673 −0,5345 −0,8018
2uut −0,2673 0,9286 −0,1429 −0,2143 .
H1 := I − =
kuk22 −0,5345 −0,1429 0,7143 −0,4286
−0,8018 −0,2143 −0,4286 0,3571
Al hacer H1 A en Octave, obtenemos la matriz
−3,7417 −1,6036
−0,0000 0,3042
A2 =
−0,0000
.
−0,3917
−0,0000 −1,0875
Notamos que, a menos de errores de redondeo (las entradas son del orden de 10−15 ), la
primera columna de esta matriz se corresponde a la de una matriz triangular superior.
A continuación, consideramos la matriz
0,3042
Ã2 := −0,3917 ;
−1,0875
(1)
en este caso, esta matriz tiene una única columna. Computamos kÃ2 k2 = 1,1952, consi-
(1)
deramos la transformación de Householder H̃2 asociada al vector u := Ã2 + kÃ2 k2 e1 , y
consideramos
1 0 ... 0
1 0 0 0
0
0 −0,2545 0,3277 0,9099
H2 := .. = .
. H̃2 0 0,3277 0,9144 −0,2377
0 0 0,9099 −0,2377 0,3401
−0,0000 −0,0000
A menos de errores de redondeo, esta es la misma matriz R que obtuvimos en el Ejemplo
5.3.1, y es un factor de la descomposición QR de A.
Por otra parte, si computamos z := H2 H1 y para el vector y de datos que definimos en el
Ejemplo 5.2.1 y corresponde al lado derecho del sistema, tenemos
−4,9443
1,1952
z= −0,4347 .
0,4313
138
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por lo tanto, para hallar la solución x ∈ R2 al problema de mı́nimos cuadrados, nos basta
con resolver el sistema triangular superior
−3,7417 −1,6036 x1 −4,9443
= .
−0,0000 −1,1952 x2 1,1952
1,7500
Es inmediato verificar que la solución es x = . Remarcamos que la obtuvimos
−1,0000
sin necesidad de computar la factorización QR de la matriz A. 4
Si bien en general no esperamos que r(x∗ ) = 0, razonemos como si esto pudiese ocurrir, y
planteemos el método de Newton-Raphson para esta función. Comenzamos de un iterado
inicial x0 ∈ Rn . Luego, dado xk ∈ Rn , linealizamos la función residuo en dicho vector,
esto es, aproximamos
r(x) ≈ r(xk ) + Jr (xk )(x − xk ).
8
Aquı́, estamos usando la expresión “no lineal” para indicar “no necesariamente lineal”. El método
que presentamos es aplicable al caso en que la dependencia respecto a algunos parámetros es lineal y
respecto a otros no lo es, como en varios de los modelos que mencionamos en el Ejemplo 5.1.1.
139
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Notemos que la función residuo r(xk ) es la resta de un término que depende de xk y otro
que no, por lo que Jr (xk ) = Jf (xk ), donde f : Rn → Rm es tal que
Para definir los iterados, igualamos la linealización del residuo al vector nulo: tomamos
xk+1 tal que
r(xk ) + Jr (xk )(xk+1 − xk ) = 0.
Como en cada paso debemos resolver un problema de mı́nimos cuadrados lineal, podemos
usar cualquiera de los métodos que analizamos en las secciones anteriores de este capı́tulo.
Al ser un método de tipo Newton-Raphson, en caso de que el problema (5.10) tenga
solución, se puede probar la convergencia si el iterado inicial está lo suficientemente cerca
de dicha solución. Asimismo, se podrı́an considerar variantes en el espı́ritu de las que
describimos en la Sección 4.8.3, como introducir amortiguamientos, evitar computar las
jacobianas en cada paso, o ir haciendo actualizaciones de rango bajo.
140
Métodos Numéricos Facultad de Ingenierı́a – IMERL
t y
0 2
1 0.7
2 0.3
3 0.1
141
Métodos Numéricos Facultad de Ingenierı́a – IMERL
142
Capı́tulo 6
Descomposición SVD
Definición 6.0.1
√ (valores singulares). Con la notación de arriba, los valores singulares
de A son σi := λi , i = 1, . . . , n. 4
143
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Teorema 6.0.1 (descomposición SVD). Dada A ∈ Mm×n (R), con rg(A) = r, existen
U ∈ Mm (R), V ∈ Mn (R) ortogonales y Σ ∈ Mm×n (R) diagonal tales que
A = U ΣV t (6.1)
y Σ es de la forma
σ1 0 . . . 0
. .
0 σ2 . . ..
Σr Or×(n−r)
Σ= , con Σr = . . ∈ Mr (R),
O(m−r)×r O(m−r)×(n−r) .. .. ... 0
0 . . . 0 σr
144
Métodos Numéricos Facultad de Ingenierı́a – IMERL
145
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Contar con la factorización SVD (6.1) de una cierta matriz A ∈ Mm×n (R) también
permite determinar de forma elegante y eficiente los espacios col(A) y ker(A). En efecto,
si rg(A) = r, entonces el espacio de columnas de A (o equivalentemente, la imagen de A)
está generado por las primeras r columnas de U ,
e incluso ese conjunto nos da una base ortonormal de col(A). Por otra parte, el núcleo de
A está dado por las últimas n − r columnas de V ,
y este conjunto es una base ortonormal de ker(A). Las funciones de Octave orth y null
usan la factorización SVD para devolver respectivamente estas bases ortonormales del
espacio de columnas y del núcleo de A.
A = U ΣV t = σ1 E1 + . . . + σn En , (6.2)
donde Ei = U (i) (V (i) )t ∈ Mm×n (R) es una matriz de rango 1. La expresión (6.2) permite
realizar aproximaciones de rango bajo de la matriz A al tomar algunos de los sumandos.
De hecho, teniendo en cuenta nuestra convención de ordenar los valores singulares en
forma decreciente, esto produce una mejor aproximación en el siguiente sentido.
Teorema
Pk 6.1.1 (Eckhart-Young). Dados A ∈ Mm×n (R) y k < r := rg(A), sea Ak :=
σ E
i=1 i i definida a partir de tomar los primeros k sumandos en (6.2). Entonces,
donde k · k2 denota la norma operador en Mm×n (R) asociada a la norma euclı́dea (Defi-
nición 2.5.1).
Este teorema tiene diversas aplicaciones. Como primer ejemplo, consideremos la compre-
sión de imágenes. Digitalmente, una imagen en blanco y negro se suele representar como
una matriz de dimensiones m × n. Cada entrada de dicha matriz representa un pı́xel,
y suele tomar valores entre 0 (negro) y 255 (blanco). La figura 6.2 muestra un ejemplo
simple de representación matricial de una imagen. El Teorema 6.1.1 nos brinda una forma
de comprimir nuestra imagen al hacer aproximaciones de rango bajo. La Figura 6.3 nos
muestra un ejemplo.
146
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 6.3: Compresión de una imagen: la original (1440 × 1920 pı́xeles) está abajo a la
derecha. Mostramos las aproximaciones dadas por el Teorema 6.1.1 de rangos k = 10
(arriba izquierda), k = 50 (arriba derecha) y k = 200 (abajo izquierda).
147
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Visualizar estos datos como puntos en R3 no necesariamente serı́a de mucha ayuda para
interpretarlos. Parece claro que las tres mediciones deben estar correlacionadas de alguna
forma: deberı́a haber una componente, llamémosla el tamaño, que permita explicar buena
parte de las diferencias entre el peso, altura y perı́metro cefálico de cada bebé. Colocamos
nuestros datos en una matriz, la centramos y hacemos que tenga varianza 1 en cada fila2
para obtener una matriz A ∈ M3×10 (R), a la que le aplicamos la factorización SVD,
A = U ΣV t . Obtenemos los valores singulares
σ1 = 4,9085, σ2 = 1,4182, σ3 = 0,9463.
Como σ1 es bastante mayor que los otros dos valores singulares, parece que buena parte
de la variabilidad de nuestros datos puede ser explicada por un solo factor (el tamaño).
Tomamos esta componente principal, es decir, la primer columna de V , y, para facilitar
la interpretación de los resultados, la normalizamos para que tenga norma infinito igual
a 1. Obtenemos el siguiente vector:
tamano = [0,90, 0,14, −0,67, 0,05, −0,09, −0,40, 1,00, −0,24, 0,27, −0,95].
2
Este paso es para evitar las distorsiones que nos podrı́a producir el hecho de que distintas mediciones
estén tomadas en unidades distintas y son de magnitudes distintas.
148
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ası́, parece ser que los bebés número 7 y 1 son los de mayor tamaño, que el último es
el que tiene menor tamaño, y si queremos ver a un bebé de tamaño “medio”, deberı́amos
buscar al número 4 o al 53 . Si queremos hacer un modelo explicativo usando solamente la
variable tamaño, consideramos la aproximación de rango 1 de nuestra matriz de datos,
E1 = σ1 U (1) (V (1) )t .
La Figura 6.4 muestra las mediciones junto a las columnas correspondientes de la matriz
E1 , que nos indica los valores predichos de acuerdo al tamaño de cada bebé. En azul se
muestran las mediciones y en rojo los valores predichos por el tamaño de cada bebé. Por
ejemplo, vemos que para nuestro modelo el bebé número 7 pesa más y mide menos de lo
que predice su tamaño, o que el bebé 9 es un poco cabezón para su tamaño.
Figura 6.4: Datos del Ejemplo 6.1.1. En azul, se representan las mediciones de cada bebé
(de arriba a abajo: peso, altura, perı́metro cefálico), y en anaranjado los valores predichos
de acuerdo al tamaño computado de cada bebé.
Además, podemos aprovechar nuestro análisis para predecir nuevos valores a partir de los
datos que obtuvimos: si conocemos una medida de un bebé, por ejemplo su peso, podemos
calcular su tamaño y luego predecir las otras mediciones. Por ejemplo, si tenemos un bebé
que pesa 10 kg, luego de restar la media de la muestra y normalizar obtenemos que tiene
un peso normalizado de −0,0401. Esta cantidad −0,0401 serı́a la primera entrada en los
datos, y como estamos haciendo una aproximación de rango 1, conocer una entrada de una
columna es suficiente para extrapolar las demás entradas. Ası́, predecimos que el tamaño
de nuestro bebé es
tamaño estimado = −0,0401/U (1, 1) = −0,0681.
3
Para nuestro conjunto de datos, el peso medio es de 10.03 kg, la altura media es de 75.8 cm, y el
perı́metro cefálico medio es de 45.7 cm.
149
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Luego, podemos “completar” los datos faltantes (en unidades normalizadas) haciendo
Vender una pelı́cula (o sistemas pay per view ). Dada una pelı́cula, determinar a qué
usuario que no la vio le podrı́a gustar.
150
Métodos Numéricos Facultad de Ingenierı́a – IMERL
este problema serı́a considerar únicamente las pelı́culas que ambos puntuaron. Como
queremos un sistema de recomendación, esto no parece una alternativa razonable.
Para remediar esto, podemos usar la regla del coseno en lugar de la distancia. Definimos
la similaridad s entre los usuarios i y j como
uti uj
s(i, j) := ,
kui k × kuj k
151
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Procedemos ahora como antes, normalizando las valoraciones de acuerdo al puntaje pro-
medio asignado por cada usuario.
152
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Notamos que con solamente usar dos valores principales, la primer fila se ya se recons-
truye razonablemente bien en los casos donde el usuario asignó valoraciones. Ahora que
tenemos nuestras predicciones, procedemos a generar recomendaciones para los usuarios.
Recomendaremos las pelı́culas que predecimos puntuarán más alto y que cada usuario no
miró.
if ~isempty(items_no_valorados)
% Ordenar los ı́tems no valorados (de mayor a menor)
valores_predichos = matriz_predicha(usuario, items_no_valorados);
[ratings_ordenados, indices] = sort(valores_predichos, ’descend’);
items_recomendados = items_no_valorados(indices);
153
Métodos Numéricos Facultad de Ingenierı́a – IMERL
end
else
fprintf(’\nEl usuario %d ya ha valorado todos los ı́tems\n’, usuario);
end
end
y ası́ sucesivamente.
De esta forma, luego de obtener la factorización SVD, el resto de los pasos se pueden
calcular a orden constante. Esta es la principal ventaja de del acercamiento a través de la
desomposición en valores singulares: escala mejor con respecto a la cantidad de datos, y
comparado al filtrado colaborativo (FC) básico que descibimos antes, resulta más eficiente
si tenemos una base de usuarios y pelı́culas grande. Resumimos las diferencias entre ambos
enfoques en la siguiente tabla.
FC Básico SVD
2
Complejidad de Preprocesamiento O(m n) O(mı́n(mn2 , m2 n))
Complejidad de Predicción O(m log m) O(k)
Complejidad de Recomendación O(mn log n) O(mn log n)
Almacenamiento Requerido O(m2 + mn) O(k(m + n))
Escalabilidad Limitada (por m2 ) Alta (con k chico)
Mejor caso de uso Pocos datos Muchos datos
Cuadro 6.1: Comparación entre el algoritmo FC básico y el algoritmo que aplica SVD.
154
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Este es un problema extremadamente simple, ya que Σ es una matriz diagonal: basta con
restringirnos a los primeros n elementos de z y dividirlos por los valores singulares de A.
Una vez que hallamos w∗ := V t x∗ , la solución de nuestro problema es x∗ = V w∗ .
La justificación de este procedimiento es análoga a la que hicimos en la Sección 5.3.1. En
efecto, usando A = U ΣV t podemos escribir la norma euclı́dea al cuadrado del residuo en
un vector x ∈ Rn como
donde en el último paso usamos que U es ortogonal. Como Σ ∈ Mm×n (R) es diagonal,
tenemos que sus últimas m − n filas son todas nulas, por lo que
Ejemplo 6.2.1 (ajuste lineal con SVD). Volvemos al Ejemplo 5.2.1/5.3.1, al que esta vez
tratamos mediante la factorización SVD de A.
5
Como rg(A) = n, sus valores singulares son todos positivos y por lo tanto Σ(1:n) es invertible.
155
Métodos Numéricos Facultad de Ingenierı́a – IMERL
>> z = U’*y;
>> w = S(1:2, 1:2) \ z(1:2);
>> x = V*w;
por lo que (At A)−1 At es una inversa a izquierda de A. Sin embargo, esta matriz en general
no es una inversa a derecha, pues
A[(At A)−1 At ] 6= Im .
Esto se deduce del hecho de que la matriz de la izquierda tiene rango menor o igual que
n, mientras que la de la derecha tiene rango igual a m > n.
156
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Para generalizar la noción de inversa de una matriz, comenzamos con el caso escalar: dado
x ∈ R \ {0}, su inverso respecto al producto es x1 . En caso de que x = 0, obviamente no
tiene inverso. Vamos a llamar pseudoinverso de x al número
+ 1/x si x 6= 0,
x :=
0 si x = 0.
A+ := V Σ+ U t .
Observación 6.2.1 (no dependencia de SVD). Se puede demostrar que la Definición 6.2.1
no depende de la factorización SVD de la matriz A. 4
Esta definición nos permite recuperar los casos conocidos en que la matriz es invertible o
en problemas de mı́nimos cuadrados con rango completo.
Demostración. Nos basta con probar la segunda afirmación, pues la primera se deduce
de ella y del hecho de que si A ∈ Mn (R) es invertible, entonces At también lo es y
(At )−1 = (A−1 )t .
6
Esta no es la definición más habitual de pseudoinversa. Se la suele definir como la única matriz
que satisface las llamadas condiciones de Moore-Penrose; referimos a [GVL13, Sección 5.5.2] para más
detalles.
157
Métodos Numéricos Facultad de Ingenierı́a – IMERL
158
Métodos Numéricos Facultad de Ingenierı́a – IMERL
x∗ = A+ y = V Σ+ U t y = V Σ+ z.
La pseudoinversa de Σ es
Σ−1
+ r Or×(m−r)
Σ := ∈ Mn×m (R),
O(n−r)×r O(n−r)×(m−r)
159
Métodos Numéricos Facultad de Ingenierı́a – IMERL
son por lo general más costosos computacionalmente que usar las ecuaciones normales,
especialmente si m n. El método clásico de ortogonalización de Gram-Schmidt que
aprendimos en Geometrı́a y Álgebra Lineal 2 no es estable computacionalmente; el método
más extendido de ortogonalización, por ser el más eficiente y preciso en general, es el de
Householder.
Finalmente, usar la factorización SVD es la opción computacionalmente más cara. Su
elevado costo viene a cambio de una gran precisión y confiabilidad en problemas delicados:
en particular, usar SVD nos permite encontrar una solución elegante en problemas de
rango deficiente.
A pesar de que no profundizamos en este punto, indicamos como referencia algunas es-
timaciones de costos computacionales de resolver un sistema con A ∈ Mm×n (R) con
distintos métodos:
mn2 n3
ecuaciones normales + método directo (Cholesky) ∼ 2
+ 3
f lops;
2n3
transformaciones de Householder + sustitución hacia atrás ∼ 2mn2 − 3
f lops;
9n3
computar SVD8 + sistema diagonal ∼ 2m2 n + 4mn2 + 2
f lops.
8
Usando el método que se describe en [QSS10, Sección 5.8.3].
160
Capı́tulo 7
Integración numérica
161
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por lo tanto, teniendo en cuenta la expresión (7.1), observamos que los pesos de la regla de
Rb
cuadratura son wi = a Lin (x) dx para todo i = 0, . . . , n: los pesos solamente dependen
de la elección de los nodos de la regla, pero no dependen de la función cuya
integral queremos aproximar. 4
A continuación, damos algunos ejemplos básicos de reglas de Newton-Cotes.
Regla del punto medio. Esta es la regla de Newton-Cotes abierta más sencilla: co-
rresponde a tomar un único punto, el punto medio del intervalo [a, b], e interpolar f por
ese punto. Obviamente, el polinomio interpolante aquı́ es la constante
a+b
p0 (x) := f .
2
162
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Regla del trapecio. En una regla de Newton-Cotes cerrada, necesariamente los dos
extremos del intervalo tienen que ser nodos. Por lo tanto, la regla de Newton-Cotes cerrada
más sencilla utiliza solamente los puntos (a, f (a)), (b, f (b)) y considera la interpolante
lineal por ellos,
f (b) − f (a)
p1 (x) := f (a) + (x − a).
b−a
La regla del trapecio consiste en la aproximación
Z b
T f (a) + f (b)
I(f ) ≈ Q (f ) := p1 (x) dx = (b − a). (7.3)
a 2
b−a
Usando la notación (7.1), tenemos x0 := a, x1 := b, y w0 = w1 := 2
.
Definición 7.1.3 (orden de una regla de cuadratura). Decimos que el orden de una regla
de integración numérica es el menor grado del polinomio que la regla no integra exacta-
mente. Esto es, una regla tiene orden m si integra exactamente a todos los polinomios de
grado menor o igual a m − 1 pero hay un polinomio de grado m tal que la regla no integra
exactamente. 4
Ejemplo 7.1.1 (órdenes de las reglas del punto medio, del trapecio, y de Simpson). Es
evidente que una regla de Newton-Cotes que use n + 1 puntos tiene que ser de orden por
lo menos igual a n + 1, ya que si f es un polinomio de grado n, al interpolarlo por n + 1
puntos obtenemos pn = f . Por lo tanto, la regla del punto medio tiene orden al menos 1,
la del trapecio tiene orden al menos 2, y la de Simpson tiene orden por lo menos 3.
163
Métodos Numéricos Facultad de Ingenierı́a – IMERL
4 4 4
3 3 3
2 2 2
1 1 1
0 0 0
0 0,5 1 1,5 2 0 0,5 1 1,5 2 0 0,5 1 1,5 2
Figura 7.1: Ilustración de las reglas del punto medio (7.2) (izquierda), del trapecio (7.3)
(centro), y de Simpson (7.4) (derecha): en ellas se aproxima la integral de la función
graficada en negro por la del polinomio interpolante graficado en rojo.
0
0 0,5 1 1,5 2
Figura 7.2: Las regla del punto medio Rtiene orden por lo menos igual a 2: si f es lineal en
b
[a, b], por simetrı́a tenemos QM (f ) = a f .
M 1 1 T f (0) + f (1) 1
Q (f ) = f = , Q (f ) = = .
2 4 2 2
164
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Una pregunta esencial que nos abocamos a responder a continuación es qué podemos decir
respecto a los errores cometidos por estas reglas de cuadratura básicas.
Teorema 7.1.1 (errores en reglas de cuadratura). Sean f : [a, b] → R integrable, I(f ) =
Rb
a
f y QM (f ), QT (f ), QS (f ) definidas por (7.2), (7.3), y (7.4), respectivamente.
a) Si f ∈ C 2 , entonces
(b − a)3 00
|QM (f ) − I(f )| ≤ kf kL∞ (a,b) , (7.5)
24
y
(b − a)3 00
|QT (f ) − I(f )| ≤ kf kL∞ (a,b) . (7.6)
12
165
Métodos Numéricos Facultad de Ingenierı́a – IMERL
b) Si f ∈ C 4 , entonces
(b − a)5 (4)
|QS (f ) − I(f )| ≤ kf kL∞ (a,b) . (7.7)
2880
f 00 (θx )
f (x) = f (xM ) + f 0 (xM )(x − xM ) + (x − xM )2 .
2
Integramos
Rb de los dos lados de esta igualdad respecto a x, y reconocemos que
a
f (x M ) dx = f (xM )(b − a) = QM (f ), para obtener
b b
f 00 (θx )
Z Z
M 0
I(f ) = Q (f ) + f (xM )(x − xM ) dx + (x − xM )2 dx.
a a 2
Rb
La primer integral del lado derecho es nula, pues a (x − xM ) dx = 0; esta es justamente
la simetrı́a que identificamos en el Ejemplo 7.1.1. Por lo tanto,
Z b 00 Z b 00
M f (θx ) 2 f (θx )
|I(f ) − Q (f )| = (x − xM ) dx ≤ (x − xM )2 dx.
a 2 a 2
Como la única información que tenemos sobre θx es que está en el intervalo (a, b), acotamos
Rb 3
|f 00 (θx )| ≤ kf 00 kL∞ (a,b) . Finalmente, es sencillo verificar que a (x − xM )2 dx = (b−a)
12
, y se
concluye (7.5).
Para demostrar (7.6), notamos que la regla del trapecio resulta de integrar la interpolante
lineal p1 por (a, f (a)), (b, f (b)), por lo que podemos usar la cota (3.13):
Es claro que las reglas de Newton-Cotes en la forma en que discutimos hasta aquı́ son de
poco uso en la práctica: si el intervalo [a, b] no es muy pequeño, entonces en general no
podemos esperar que las reglas del punto medio, del trapecio, o de Simpson produzcan una
166
Métodos Numéricos Facultad de Ingenierı́a – IMERL
buena aproximación a la integral de f en [a, b]. Esto es análogo a usar una interpolante de
grado bajo en un intervalo grande: seguramente dé una aproximación de poca calidad a la
función que queremos interpolar. Una opción que puede resultar natural es aumentar el
grado polinomial; sin embargo, si recordamos nuestra discusión en la Sección 3.3.2 sobre
el fenómeno de Runge, rápidamente caemos en la cuenta de que hacer interpolaciones de
Newton-Cotes –que usan nodos equiespaciados– no es una buena idea a menos que las
funciones cuyas integrales queremos estimar sean “buenas”, como la del Ejemplo 3.3.1.
Otra propiedad que vuelve no deseable a las reglas de Newton-Cotes de orden alto es
que en ellas suele haber pesos con signo negativo. Esto puede conducir a sumas largas en
las que hay varios términos de signos opuestos y, en consecuencia, las hacen pasibles de
cancelaciones numéricas.
Podemos considerar las mismas alternativas que en la Sección 3.3.2: o elegimos mejor
los nodos de cuadratura o nos mantenemos con una interpolación de grado bajo pero
dividimos [a, b] en subintervalos pequeños. Respecto a la primer alternativa, una opción
relacionada con los nodos de Chebyschev es la llamada cuadratura de Clenshaw-Curtis,
sobre la que no vamos a profundizar; otra opción, que exploramos en la Sección 7.4, es la
llamada cuadratura gaussiana. La segunda alternativa da lugar a las reglas compuestas,
que describimos a continuación.
hi := xi+1 − xi
y, para ∗ ∈ {M, T, S}1 , ponemos Q∗i (f ) al valor resultante de aplicar la regla de cuadratura
∗ en el intervalo [xi , xi+1 ]. Como las integrales son aditivas respecto a los intervalos,
tenemos
Z b n−1 Z xi+1
X n−1
X
I(f ) = f (x) dx = f (x) dx ≈ Q∗i (f ) =: Q∗,c (f ).
a i=0 xi i=0
1
Esta notación solamente quiere indicar que vamos a considerar las reglas del punto medio, del trapecio,
y de Simpson.
167
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Si usamos la notación (7.1), tenemos que los nodos de la regla son x0 , . . . , xn , y los pesos
son w0 = h20 , wi = hi−12+hi para i = 1, . . . , n − 1, y wn = hn−1
2
.
Observación 7.2.1 (notación). Nuestra notación para integrales compuestas no hace re-
ferencia a la partición que hayamos tomado, pero obviamente el valor que calculemos
depende de ella. 4
Teorema 7.2.1 (errores de reglas compuestas). Sean f : [a, b] → R integrable, a = x0 <
Rb
. . . < xn = b, hi = xi+1 − xi (i = 0, . . . , n − 1). Llamamos I(f ) = a f y QM,c (f ),
QT,c (f ), QS,c (f ) a las reglas del punto medio, del trapecio, y de Simpson compuestas,
respectivamente.
a) Si f ∈ C 2 , entonces
n−1 3
X h
|Q M,c
(f ) − I(f )| ≤ i
kf 00 kL∞ (xi ,xi+1 ) ,
i=0
24
y
n−1 3
X h
|Q T,c
(f ) − I(f )| ≤ i
kf 00 kL∞ (xi ,xi+1 ) . (7.9)
i=0
12
b) Si f ∈ C 4 , entonces
n−1
X h5i
|QS,c (f ) − I(f )| ≤ kf 00 kL∞ (xi ,xi+1 ) .
i=0
2880
En caso de que las particiones sean uniformes, podemos escribir una cota de error en
forma más compacta.
Corolario 7.2.2 (particiones uniformes). En las mismas hipótesis que en el Teorema
7.2.1, si asumimos además que hi = h := (b−a)
n
para todo i = 0, . . . , n, tenemos:
168
Métodos Numéricos Facultad de Ingenierı́a – IMERL
a) Si f ∈ C 2 , entonces
(b − a)h2 00
|QM,c (f ) − I(f )| ≤ kf kL∞ (a,b) ,
24
y
(b − a)h2 00
|QT,c (f ) − I(f )| ≤ kf kL∞ (a,b) .
12
b) Si f ∈ C 4 , entonces
(b − a)h4 (4)
|QS,c (f ) − I(f )| ≤ kf kL∞ (a,b) .
2880
Demostración. Las desigualdades se deducen inmediatamente del teorema anterior, usan-
do que para todo i = 0, . . . , n − 1 y toda función g : [a, b] → R, se cumple kgkL∞ (xi ,xi+1 ) ≤
kgkL∞ (a,b) , y que
n−1
(b − a) X b−a
h= ⇒ 1=n= .
n i=0
h
2(2/n)2 00 100
|I(f ) − QT,c (f )| ≤ kf kL∞ (a,b) = 2 .
12 3n
Por ejemplo, si queremos garantizar que el error en nuestra aproximación sea menor que
10−4 , basta con tomar n tal que se cumpla
100 103
< 10−4 ⇒ n > √ ' 577,35,
3n2 3
169
Métodos Numéricos Facultad de Ingenierı́a – IMERL
por lo que alcanza con tomar n = 578 subintervalos. Al implementar la regla del trapecio
compuesta en Octave, nos encontramos con que el error en ese caso es aproximadamente
1,48 × 10−7 . Si bien esto está de acuerdo con la teorı́a, ya que demostramos que tomar 578
subintervalos es suficiente para que el error de integración sea menor que 10−4 , da la sensa-
ción de que nuestra cota no es muy precisa. Tomando 578 subintervalos estamos refinando
mucho más de lo necesario. De hecho, experimentando en Octave, uno se encuentra con
que basta con solamente tomar n = 14 para que el error de integración satisfaga nuestro
requerimiento de precisión: en ese caso, el error es aproximadamente 6,06 × 10−5 . 4
7.3. Adaptatividad
Sigamos interpretando los resultados del Ejemplo 7.2.1. Analizando un poco más en detalle
cómo llegamos a las estimaciones del Corolario 7.2.2, observamos que acotamos la norma
L∞ de derivadas de f (en el ejemplo, las derivadas segundas) en cada subintervalo por la
norma L∞ de dicha derivada en todo el intervalo [a, b]. La Figura 7.3 muestra el gráfico
de la derivada segunda de la función de Runge en el intervalo [0, 1] (notar que f 00 es una
función par): la cota |f 00 (x)| ≤ 50 es relevante solamente para x cercano al origen. Ya
cuando x ≈ 1/2, tenemos |f 00 (x)| ≈ 2. Tomar subintervalos tan pequeños para x ≥ 1/2 es
demasiado para el requerimiento de precisión que tenemos. En este ejemplo es claro porque
conocemos y podemos graficar la derivada segunda de f y sabemos que f es “mucho más
lineal” cuando x está alejado de 0 que cuando está cerca.
En esta sección, mostramos algunas ideas sobre adaptatividad: nuestro objetivo es ha-
cer reglas de cuadratura compuestas que de alguna forma puedan determinar cuándo es
170
Métodos Numéricos Facultad de Ingenierı́a – IMERL
necesario tomar intervalos pequeños y cuándo no vale la pena hacerlo. El objetivo es,
dada una cierta tolerancia, computar la integral de f sobre el intervalo [a, b] con esa to-
lerancia y de forma eficiente. Para ilustrar las ideas principales, vamos a fijar una regla
de cuadratura: solamente vamos a considerar que nuestra aproximación se hace mediante
la regla del trapecio compuesta. El desafı́o está en encontrar una partición del intervalo
[a, b] que nos permita alcanzar nuestro objetivo.
Supongamos que implementamos una regla del trapecio compuesta sobre una partición
a = x0 < . . . < xn = b. Descomponemos el error de integración numérica como la suma
de los errores locales, asociados a la integración en cada subintervalo:
n−1
X
T,c
Q (f ) − I(f ) = QTi (f ) − Ii (f ),
i=0
Rx
donde QTi (f ) está dado por (7.8) e Ii (f ) := xii+1 f . Sea ε nuestra tolerancia: queremos
que sea |QT,c (f ) − I(f )| ≤ ε, e idealmente no mucho menor, ya que si lo fuera querrı́a
decir que seguramente nuestro algoritmo refinó al intervalo más de lo necesario.
En principio, no sabemos cómo ni en cuántas partes tenemos que partir el intervalo [a, b]
para que el error de integración sea menor o igual a ε. Para empezar, supongamos que
podemos calcular los errores locales exactamente. Si comenzamos nuestro algoritmo con
una regla del trapecio en todo el intervalo [a, b], entonces podemos comparar si el error es
menor a ε o no. Si lo es, terminamos; si no lo es, tenemos que seguir partiendo nuestro
intervalo. Si partimos el intervalo [a, b] a la mitad, ahora tendremos dos subintervalos:
¿cómo deberı́a ser el error local en cada uno de esos subintervalos para que el algoritmo
termine? Una respuesta es que en cada uno de esos subintervalos el valor absoluto del
error local sea menor que ε/2: como tenemos dos subintervalos,
|QT,c (f ) − I(f )| = |(QT0 (f ) − I0 (f )) + (QT1 (f ) − I1 (f ))|
≤ |QT0 (f ) − I0 (f )| + |QT1 (f ) − I1 (f )| < ε.
El factor 1/2 multiplicando a ε en esta tolerancia local está relacionado con el hecho de
que los dos subintervalos miden la mitad que el intervalo [a, b]. En efecto, si logramos que
el error local en un subintervalo sea como ε multiplicado por la proporción del intervalo
[a, b] que el subintervalo cubre, entonces habremos terminado: si para cada i = 0, . . . , n−1,
se tiene
hi
|QTi (f ) − Ii (f )| ≈ ε, (7.10)
b−a
entonces
n−1 n−1 n−1
X X ε X
|QT,c (f ) − I(f )| = QTi (f ) − Ii (f ) ≤ |QTi (f ) − Ii (f )| ≈ hi = ε,
i=0 i=0
b − a i=0
Pn−1
porque i=0 hi = b − a. La estimación de arriba es lo mejor que podemos lograr si no
tenemos un control sobre los signos de los errores, lo que es bastante difı́cil de lograr en
la práctica.
171
Métodos Numéricos Facultad de Ingenierı́a – IMERL
De este modo, comenzamos a delinear cómo serı́a una cuadratura adaptativa usando la
regla del trapecio compuesta: dadas una función f : [a, b] → R integrable una tolerancia
ε > 0, y un intervalo I ⊂ [a, b] de longitud `,
172
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ejemplo 7.3.1. Implementamos nuestro Algoritmo adaptativo usando la regla del tra-
pecio compuesta a la función que consideramos en la Sección 1.2.1: queremos computar
Z 1
I= esen x dx.
0
Ponemos una tolerancia ε = 10−4 . La Figura 7.4 (izquierda) muestra el gráfico del inte-
grando f (x) = esen x en el intervalo y los puntos en los que el algoritmo hizo subdivisiones.
En total, se utilizaron 57 nodos, y obtuvimos el valor
I ≈ 1,6319.
Este valor es consistente con lo que estimamos en el Ejemplo 1.2.1 usando sumas de
R1
Figura 7.4: Cuadratura adaptativa. Izquierda: cálculo de 0 esen x dx con tolerancia ε =
R5 3x
10−4 ; derecha: cálculo de 0 2x−sen
1+x5
dx con tolerancia ε = 10−2 .
Riemann, pero con bastante más particiones del intervalo. En este ejemplo, notamos que
el gráfico de f (x) parece ser “bastante lineal”2 por lo que los subintervalos tomados son
bastante uniformes.
El efecto adaptativo de nuestro algoritmo es más evidente cuando consideramos integran-
dos cuya derivada segunda
R 5 varı́a mucho. La Figura 7.4 (derecha) muestra los intervalos
3x
tomados para calcular 0 2x−sen
1+x5
dx con una tolerancia ε = 10−2 . La partición es mucho
más fina cuando x está aproximadamente entre 1/2 y 1 que para x cercano a 5. 4
173
Métodos Numéricos Facultad de Ingenierı́a – IMERL
y buscamos w1 , w2 ∈ R, x1 , x2 ∈ [−1, 1] tales que la regla sea exacta para todos los
polinomios de grado menor o igual que 3. Para este fin, basta con elegir nuestra base
favorita del espacio de polinomios de grado menor o igual que 3, imponer que la regla sea
exacta en cada polinomio de dicha base, y resolver el problema resultante. Tomamos la
base monomial {1, x, x2 , x3 }, y queremos que se cumplan:
Z 1
2= 1 dx = w1 1 + w2 1 ⇒ w1 + w2 = 2,
−1
Z 1
0= x dx = w1 x1 + w2 x2 ⇒ w1 x1 + w2 x2 = 0,
−1
Z 1 (7.11)
2/3 = x2 dx = w1 x21 + w2 x22 ⇒ w1 x21 + w2 x22 = 2/3,
−1
Z 1
0= x3 dx = w1 x31 + w2 x32 ⇒ w1 x31 + w2 x32 = 0.
−1
Tenemos que resolver el sistema (7.11). Multiplicamos la segunda ecuación por x22 y se la
restamos a la última para obtener
174
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Observación 7.4.1. El Ejemplo 7.4.1 es un caso pequeño y que podemos tratar a mano.
Para problemas más grandes, podemos utilizar métodos iterativos para resolver numéri-
camente los sistemas de ecuaciones resultantes, como los que discutimos en la Sección 4.8.
La cuadratura gaussiana está relacionada con la teorı́a de polinomios ortogonales y exis-
ten algoritmos que explotan este vı́nculo para computar pesos y nodos en forma eficiente
y estable. Al lector interesado en el tema, le referimos a [QSS10, Capı́tulo 10] para más
detalles. 4
Observación 7.4.2 (otros intervalos). Las cuadraturas gaussianas se suelen expresar sobre
el intervalo [−1, 1]. Cualquier integral sobre un intervalo [a, b] se la puede llevar a una
integral sobre el intervalo [−1, 1] mediante un cambio de variables afı́n, y los pesos y
nodos cambian correspondientemente, como se discutirá en el práctico. 4
Observación 7.4.3 (pesos son positivos). Desde el punto de vista computacional, una
propiedad importante de la cuadratura gaussiana es que, en general, los pesos que se
obtienen son positivos. Se puede encontrar una demostración en [SM03, Sección 10.2]. 4
3
Formalmente, también existe la solución x1 = √13 , w1 = 1, x2 = − √13 , w2 = 1, pero es claro que
esta solución da lugar a la misma regla de cuadratura que la que escribimos arriba.
175
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Si bien los términos del lado derecho en (7.14) son iguales a los de (7.11), la resolución de
este sistema parece más compleja, y no parece claro que los nodos y pesos se deban dis-
tribuir simétricamente como en el Ejemplo 7.4.1. Implementamos en Octave el método de
Newton-Raphson que describimos en la Sección 4.8.2, y obtenemos los valores numéricos
x̂1 ' 0,2476, x̂2 ' 0,8192, ŵ1 ' 0,7131, ŵ2 ' 1,0051.
176
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Usamos el sı́mbolo de tilde ˆ para denotar que estos nodos y pesos son los correspondientes
al intervalo [0, 1]. Si queremos una regla para (7.13) en el intervalo [0, a], simplemente
hacemos el cambio de variables x̂ = x/a y aplicamos la regla de cuadratura hallada a la
función fˆ(x̂) := f (ax̂),
Z a Z 1 Z 1
x/a
f (x) e dx = a x̂
f (ax̂) e dx = a fˆ(x̂) ex̂ dx
0
0 0
ˆ ˆ
' a ŵ1 f (x̂1 ) + ŵ2 f (x̂2 ) = aŵ1 f (ax̂1 ) + aŵ2 f (ax̂2 ).
177
Métodos Numéricos Facultad de Ingenierı́a – IMERL
178
Capı́tulo 8
179
Métodos Numéricos Facultad de Ingenierı́a – IMERL
entonces podemos definir las n variables auxiliares y1 (t) = y(t), y2 (t) = y 0 (t) = y10 (t), . . . ,
0
yn (t) = yn−1 (t), y la ecuación (8.3) puede escribirse como una ecuación de primer orden
para la función vectorial y : R → Rn , de la forma
y10 (t) y2 (t)
y0 (t) y3 (t)
2
y0 (t) = .. = =: f (t, y(t)).
..
. .
0
yn (t) f (t, y1 (t), y2 (t), . . . , y(n−1) (t))
En el Ejemplo 8.1.3 aplicamos esta idea a una ecuación diferencial de segundo orden
proveniente de aplicar las leyes de Newton.
A continuación mostramos algunos ejemplos relevantes de ecuaciones diferenciales.
Ejemplo 8.1.1 (problema test). Consideremos f (t, y(t)) = λy(t) en (8.2), donde λ ∈ R
es un dato. Nuestra ecuación diferencial tiene entonces la forma y 0 (t) = λy(t). Es sencillo
verificar que, para todo C ∈ R, la función
y(t) = Ceλt
180
Métodos Numéricos Facultad de Ingenierı́a – IMERL
las derivadas segundas de x. Sin embargo, podemos escribirla en la forma (8.2) de un modo
sencillo. Consideremos la velocidad horizontal de la masa, y(t) := x0 (t). Como y 0 (t) =
x00 (t), podemos reescribir (8.4) usando ambas variables (desplazamiento y velocidad) como
0
x (t) = y(t),
y 0 (t) = − m
k
x(t).
2 x(t)
Luego, podemos definir una función incógnita vectorial y : R → R , y(t) := y
y(t)
escribir nuestra ecuación diferencial como
0 0 1
y (t) = Ay(t), donde A = k .
−m 0
4
2
Esta es la llamada ley de Hooke.
181
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Ejemplo 8.1.4 (problema de valores iniciales). Volvemos al Ejemplo 8.1.1. Sabemos que
las funciones que cumplen la ecuación diferencial y 0 (t) = ay(t) son de la forma
y encontramos la solución
y(t) = y0 ea(t−t0 ) .
4
182
Métodos Numéricos Facultad de Ingenierı́a – IMERL
entonces existe > 0 tal que el problema (8.6) tiene una solución única en I := [t0 , t0 + ).
Ejemplo 8.1.5 (necesidad de la hipótesis de ser Lipschitziana). Consideremos el problema
de valores iniciales ( p
y 0 (t) = |y(t)|,
y(0) = 0.
Es sencillo verificar que las funciones
0 si t ≤ 0
y1 (t) ≡ 0 e y2 (t) = t2
si t > 0
4
son soluciones del problema,
p y por lo tanto no tenemos unicidad. Esto se justifica porque
la función f (t, y) = |y| no está en las hipótesis del Teorema de Picard. 4
En lo que sigue, vamos a asumir que los problemas que estamos tratando están en las
hipótesis del Teorema de Picard. Vamos a asumir que nuestros problemas tienen solución
única y que ésta está definida en algún intervalo I conocido.
En la teorı́a de las ecuaciones diferenciales, se suele usar la palabra estabilidad para indicar
un concepto para el que en nuestro curso hemos empleado el término condicionamiento:
informalmente, esta noción corresponde a que pequeños cambios en el valor inicial y0
induzcan pequeños cambios en los valores de y(t) para todo t ≥ t0 .
Definición 8.1.2 (estabilidad de una solución de un problema de valores iniciales). La
solución y(t) del problema de valores iniciales (8.6) se dice estable (en el futuro) si para
todo > 0 existe δ > 0 tal que si ŷ(t) satisface la misma ecuación diferencial pero con
condición inicial ŷ(t0 ) = ŷ0 con kŷ0 − y0 k ≤ δ, entonces
183
Métodos Numéricos Facultad de Ingenierı́a – IMERL
nuestra experiencia en este curso y en particular del Capı́tulo 1, sabemos que los errores
de redondeo son una parte de la historia, y que otra parte fundamental son los errores de
truncamiento, asociados a la aproximación de la ecuación diferencial mediante métodos
numéricos. En la Sección 8.3 damos nociones de estabilidad de métodos numéricos para
la resolución de problemas de valores iniciales, que contemplan el comportamiento de los
errores de truncamiento.
Comencemos nuestro análisis de métodos numéricos para la aproximación de soluciones.
Consideremos un problema de valores iniciales de la forma (8.6),
(
y0 (t) = f (t, y(t))
.
y(t0 ) = y0
Nuestro objetivo es encontrar una secuencia {(tk , yk )}k≥1 (que puede ser finita o infinita,
dependiendo de si en (8.6) es I = [t0 , Tmax ) o I = [t0 , ∞)) de forma tal que yk ≈ y(tk )
para todo k. Los pasos hk := tk+1 − tk pueden ser elegidos a priori o ser determinados a
partir de ciertos requisitos de tolerancia en la solución computada.
El teorema fundamental del cálculo muestra que, si y es solución al problema de valores
iniciales (8.6), entonces debe satisfacer
Z t+h
y(t + h) = y(t) + f (s, y(s)) ds.
t
Para el caso especial en el que f no depende de y, la expresión (8.1.1) nos permite escribir
Z tn+1
yn+1 = yn + f (s) ds,
tn
y por lo tanto, obtenemos los valores de {yk } resolviendo por cuadraturas mediante cual-
quiera de los métodos vistos en el Capı́tulo 7. Sin embargo, no podemos usar cuadraturas
directamente para resolver este problema si la función f depende de y, dado que no cono-
cemos la función y(s) y, por lo tanto, no podemos evaluar el integrando. Para determinar
{yk } en el caso general, se reemplaza la ecuación diferencial por una o varias ecuaciones
algebraicas. La diferencia entre los diferentes métodos de resolución de ecuaciones ordina-
rias radica precisamente en qué forma toman las ecuaciones algebraicas para determinar
los iterados.
En la Sección 8.2 comenzamos el estudio teórico de los métodos numéricos para problemas
de valores iniciales con los métodos más básicos, los llamados métodos de Euler. La Sección
8.3 introduce los conceptos fundamentales para el análisis de métodos numéricos: consis-
tencia, estabilidad y convergencia. Aplicamos estos conceptos a un nuevo método, el del
trapecio, en la Sección 8.4. Posteriormente, en la Sección 8.5 consideramos y mostramos
cómo se construye una familia de métodos, los llamados de Runge-Kutta; estos métodos
tienen ciertas caracterı́sticas que hacen que sean la base de un buen número de los solvers
de Matlab/Octave, y explicamos las ideas clave del funcionamiento de uno de ellos. La
184
Métodos Numéricos Facultad de Ingenierı́a – IMERL
185
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Luego, para resolver el problema (8.6) numéricamente con el método de Euler hacia atrás,
tomamos y0 ∈ Rn y luego realizamos la iteración
Observación 8.2.1 (una diferencia fundamental). Si bien los métodos (8.9) y (8.10) son en
apariencia muy similares, tienen diferencias fundamentales. A primera vista, es claro que si
conocemos yk entonces para determinar el iterado siguiente yk+1 con el método de Euler
hacia adelante, basta con reemplazar los valores conocidos de tk , yk en el lado derecho
de (8.9). Esto justifica la denominación de explı́cito para ese método. En cambio, para
determinar yk+1 con el método de Euler hacia atrás, nos encontramos con una ecuación
que debemos resolver. Dependiendo de cómo sea la función f , es posible que debamos usar
alguno de los métodos para ecuaciones no lineales que tratamos en el Capı́tulo 4. Ası́, en
el método de Euler hacia atrás en general el iterado yk+1 queda definido implı́citamente.
En un método numérico para discretizar ecuaciones diferenciales ordinarias, la diferencia
entre ser explı́cito o implı́cito tiene consecuencias tanto prácticas como analı́ticas. Desde el
punto de vista práctico, en general es mucho más sencillo computar iterados con métodos
explı́citos, y el costo computacional de agregar una resolución de una ecuación no lineal en
cada paso puede ser prohibitiva en varias aplicaciones. Sin embargo, los métodos implı́citos
por lo general gozan de mejores propiedades de estabilidad que los explı́citos, lo que a su
vez permite tomar pasos más largos y, en consecuencia, tomar menos pasos. Profundizamos
sobre este último punto en la Sección 8.3.3. 4
Notamos que, en esta iteración, podemos tomar los pasos sin mayor dificultad: dado yk ,
basta con evaluar el producto que tenemos del lado derecho para determinar yk+1 .
186
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Por otra parte, si aplicamos el método de Euler hacia atrás (8.10) para este problema,
obtenemos
yk+1 = yk + ayk+1 (1 − yk+1 )hk . (8.12)
Aquı́ observamos que no podemos evaluar el lado derecho directamente, sino que yk+1 está
dado implı́citamente. En este caso, como obtenemos una función cuadrática, podemos
despejar
p
2 −(1 − ah k ) ± (1 − ahk )2 + 4ahk yk
ahk yk+1 + (1 − ahk )yk+1 − yk = 0 ⇒ yk+1 = .
2ahk
Tı́picamente se hacen solamente unas pocas iteraciones en j antes de dar el siguiente paso
∂f ∂f
en k. Observamos que Jg (y) = hk ∂y (tk+1 , y), donde ∂y denota la matriz de derivadas
parciales de f respecto a la variable y. Por lo tanto, si hk es lo suficientemente pequeño,
187
Métodos Numéricos Facultad de Ingenierı́a – IMERL
podemos garantizar que ρ(Jg (y)) < 1 y tendremos que la iteración es localmente con-
vergente. Además, a menos que y varı́e dramáticamente rápido, el iterado inicial yk va a
estar cerca de yk+1 .
Por ejemplo, apliquemos esta iteración a (8.12) para tomar el primer paso con los valores
de parámetros a = 1, y0 = 0,1, t0 = 0, y paso constante h = 0,1. Nuestra iteración toma
y10 = 0,1, e
y1j+1 = y0 + ay1j (1 − y1j )h1 = 0,1 + 0,1y1j (1 − y1j ), j ≥ 0.
Esto da lugar a los iterados
188
Métodos Numéricos Facultad de Ingenierı́a – IMERL
aquı́ consideramos métodos que para generar yk+1 solamente hacen uso del valor de yk ,
en principio es posible también utilizar la información respecto a iterados anteriores. Esto
conduce a la primera definición.
Observación 8.3.1 (métodos de Euler). Los métodos de Euler (8.9) y (8.10) son ambos de
paso simple. El método de Euler hacia adelante es explı́cito y el de Euler hacia atrás es
implı́cito, tal como adelantamos en la Observación 8.2.1. 4
Este es un método multipaso porque para determinar yk+1 requerimos los valores de yk
e yk−1 , y es un método explı́cito porque para calcular yk+1 basta con reemplazar con los
valores de yk e yk−1 y operar.
8.3.2. Consistencia
Al analizar qué tan cerca está la solución que obtenemos mediante un método numérico
de la solución exacta de un problema de valores iniciales, hay dos nociones de error rele-
vantes. Una ecuación diferencial ordinaria admite una familia de soluciones y, al fijar una
condición inicial –esto es, al considerar un problema de valores iniciales–, seleccionamos
una de ellas. En nuestro método numérico para un problema de valores iniciales, al ir
haciendo pasos, a menos que el valor que encontremos sea exacto nos iremos cambiando
de solución dentro de esa familia. La Figura 8.3 ilustra este hecho: comenzando de (t0 , y0 ),
que se encuentra sobre la curva solución, al tomar el primer paso obtenemos un punto
(t1 , y1 ) que está sobre una curva distinta de la familia de soluciones. Al tomar otro paso,
189
Métodos Numéricos Facultad de Ingenierı́a – IMERL
nos vamos a mover a otra curva distinta, y ası́ sucesivamente. En el instante t3 , nuestro
método computa un valor y3 . Nos interesa conocer la discrepancia entre y3 e y(t3 ), la
distancia entre los puntos en rojo y verde hacia la derecha de la figura. Como el análisis
de los métodos se suele realizar considerando lo que ocurre al tomar solamente un paso, va
a ser útil considerar la discrepancia entre y3 y la solución exacta al problema que tiene la
misma ecuación diferencial pero condición inicial y(t2 ) = y2 en el instante t3 : tal solución
exacta es lo que llamamos u(t3 ) y marcamos con un punto en azul en la figura.
y(t3 )
u(t3 )
y3 y3
y2
y1
y0
t0 t1 t2 t3
Definición 8.3.3 (errores global y local). Sea {(tk , yk )} dado por un método numérico
para aproximar el problema (8.6), cuya solución denotamos por y. Los errores globales
son {ek } dados por
ek := yk − y(tk ). (8.15)
Por otra parte, los errores locales de truncamiento (o errores de consistencia) {``k }
se definen de la siguiente manera: ` 0 = 0 y, para todo k ≥ 0,
4
Definición 8.3.4 (consistencia). Diremos que un método numérico es consistente si
para todo k ≥ 0 se cumple
`k+1
→ 0 con hk = tk+1 − tk → 0.
hk
190
Métodos Numéricos Facultad de Ingenierı́a – IMERL
La elección de qué norma usar aquı́ es irrelevante porque todas las normas en Rn son
equivalentes. 4
Observación 8.3.2 (no es un error de tipeo). Puede parecer llamativo que, en la definición
anterior, para referirnos a un método de orden p, en la fórmula (8.17) aparece un exponente
p + 1. La razón es que para definir el orden de un método usamos los errores locales de
truncamiento. La cantidad que nos interesa es en realidad el error global y, en el pasaje de
local a global tendremos un efecto de acumulación de error que, en el mejor de los casos,
nos hace perder un orden (ver (8.20) más abajo). 4
Ejemplo 8.3.1 (los métodos de Euler son de primer orden). Analicemos los errores de
consistencia para los métodos de Euler explı́cito e implı́cito. Para ello, dado yk , debemos
estudiar los errores de truncamiento dados por (8.16).
Para el método de Euler hacia adelante (8.9) tenemos
8.3.3. Cero-estabilidad
La noción de estabilidad de un método numérico para resolver una ecuación diferencial
ordinaria es análoga a la de estabilidad de soluciones: queremos que pequeñas pertur-
baciones en un determinado instante no generen que la solución que obtengamos crezca
indefinidamente en tiempos posteriores.
El concepto de cero-estabilidad garantiza que en un intervalo acotado fijo pequeñas pertur-
baciones de los datos producen perturbaciones acotadas de la solución numérica, cuando
la longitud del paso máximo tiende a 0. El requerimiento de estabilidad de un método
191
Métodos Numéricos Facultad de Ingenierı́a – IMERL
aparece, en primer lugar, de la necesidad de mantener bajo control los errores introducidos
al estar utilizando aritmética de punto flotante. En efecto, si un método numérico no es
cero-estable, los errores de redondeo tanto al tomar y0 como en el proceso de evaluar f
harı́an que la solución computada fuera completamente inútil.
kyk − zk k ≤ Cky0 − z0 k ∀k = 0, . . . , K.
8.3.4. Convergencia
La pregunta de convergencia refiere a qué ocurre cuando en nuestro método numérico
hacemos que la longitud de todos los pasos tienda a 0. Tı́picamente, en la práctica, incluso
cuando el problema (8.6) esté definido en un intervalo no acotado I = [t0 , ∞), al resolverlo
numéricamente y hacer que la longitud máxima del paso tienda a 0 nos interesa controlar
los errores hasta algún tiempo máximo finito. Por lo tanto, en esta sección vamos a
trabajar con problemas sobre intervalos finitos I = [t0 , Tmax ).
Definición 8.3.6 (convergencia). Diremos que un método numérico para resolver el pro-
blema (8.6) es convergente si, para todo T ∈ (t0 , Tmax ) fijo, si tomamos t0 < t1 < . . . <
tK = T arbitrarios y llamamos h := máxk=0,...,K−1 hk = máxk=0,...,K−1 tk+1 − tk , tenemos
En palabras, la convergencia de un método quiere decir que, para todo tiempo fijo T , al
hacer que los largos del paso tiendan a 0, se tiene que el valor computado por el método
en tiempo T converge al valor de la solución en dicho tiempo. 4
192
Métodos Numéricos Facultad de Ingenierı́a – IMERL
El teorema anterior tiene un significado preciso en el contexto que hemos discutido hasta
ahora, pero su significado “filosófico” también es válido en otros contextos. Muchas veces,
la convergencia de métodos numéricos (que el lı́mite de lo que computamos sea la función
buscada) es equivalente a alguna forma de consistencia (que las ecuaciones que se quiere
aproximar se verifiquen en el lı́mite en que los parámetros de discretización tienden a
cero) más alguna forma de estabilidad (que el método no magnifique descontroladamente
los errores de redondeo y de aproximación numérica).
Recordamos que los métodos de un paso consistentes son cero-estables, por lo que el
Teorema 8.3.2 establece su convergencia. Podemos hacer más preciso el contenido de ese
teorema en ese caso, vinculando los errores de truncamiento y globales.
193
Métodos Numéricos Facultad de Ingenierı́a – IMERL
La solución de este problema es y(t) = eλt , que cumple y(t) → 0 con t → ∞ porque
λ < 0. Consideremos un método numérico con paso fijo h > 0, que produce una sucesión
{(tk , yk )}. Deseamos que yk → 0 con k → ∞: si esto ocurre, decimos que el método es
absolutamente estable3 para ese paso.
Analicemos la estabilidad de los métodos de Euler, comenzando por el explı́cito. Si apli-
camos (8.9) a (8.21), nos encontramos con que, para cada k ≥ 0,
194
Métodos Numéricos Facultad de Ingenierı́a – IMERL
195
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Para que w → 0 con t → ∞, necesitamos que el radio espectral de esta matriz jacobiana
sea menor que 1. Como los valores propios de una matriz real en general son complejos,
es matemáticamente más correcto permitir λ ∈ C en el problema test para definir una
región de estabilidad en lugar de un intervalo.
Formalmente, para el problema test y 0 = eλy , la región de estabilidad es el conjunto
A := {z := hλ ∈ C : |yk | → 0}.
196
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 8.4: Solución a (8.21) para λ = −2 en I = [0, 20]. En rojo se representa la solución
exacta, mientras que en azul se muestran las soluciones producidas por los métodos de
Euler hacia adelante (primera fila) y de Euler hacia atrás (segunda fila). En todos los
casos se tomaron pasos uniformes, con valores h = 2/5, h = 2/3, h = 1 y h = 4/3, de
izquierda a derecha.
4
En el contexto de ecuaciones en derivadas parciales, este método también recibe el nombre de método
de Crank-Nicolson.
197
Métodos Numéricos Facultad de Ingenierı́a – IMERL
resulta al combinar los métodos de Euler. En efecto, para el problema de valores iniciales
(8.6), este método consiste en aproximar y(tk+1 ) por
f (tk , yk ) + f (tk+1 , yk+1 )
yk+1 := yk + hk . (8.23)
2
En otras palabras, los iterados del método del trapecio se calculan poniendo del lado
derecho la mitad de la suma entre los lados derechos correspondientes a los iterados de
los métodos de Euler hacia adelante y de Euler hacia atrás.
A partir de la definición (8.23), se infiere inmediatamente que el método del trapecio es
implı́cito.
Proposición 8.4.1. El método del trapecio es incondicionalmente absolutamente estable
y es de segundo orden.
Queremos probar que el error de truncamiento cumple ` k+1 = yk+1 − u(tk+1 ) = O(h3k ),
asumiendo que u es de clase C 3 . Para ello, hacemos un desarrollo de Taylor alrededor de
tk , evaluado en tk+1 :
h2k 00
u(tk+1 ) = u(tk ) + hk u0 (tk ) + u (tk ) + O(h3k )
2 (8.24)
h2k 00 3
= yk + hk f (tk , yk ) + u (tk ) + O(hk ).
2
00
Para tratar el término u (tk ), recordamos de la Sección 1.5.3 que la aproximación por
diferencias hacia adelante es de primer orden, esto es,
00 u0 (tk+1 ) − u0 (tk )
u (tk ) = + O(hk ).
hk
198
Métodos Numéricos Facultad de Ingenierı́a – IMERL
h2k 00
u(tk+1 ) = yk + hk f (tk , yk ) + u (tk ) + O(h3k )
2
hk
= yk + hk f (tk , yk ) + [f ((tk+1 , u(tk+1 )) − f (tk , yk )] + O(h3k )
2
hk
= yk + [f ((tk+1 , u(tk+1 )) + f (tk , yk )] + O(h3k ).
2
Teniendo en cuenta la definición (8.23) de los iterados con el método del trapecio, obte-
nemos
hk
` k+1 = [f (tk+1 , yk+1 ) − f (tk+1 , u(tk+1 ))] + O(h3k ). (8.25)
2
Como f es Lipschitziana respecto a la segunda variable, existe una constante C > 0 tal
que
kf (tk+1 , yk+1 ) − f ((tk+1 , u(tk+1 ))k ≤ Ckyk+1 − u(tk+1 )k = Ck``k+1 k.
Por lo tanto, tomando normas en (8.25) y usando la desigualdad triangular, deducimos
hk
k``k+1 k ≤ C k``k+1 k + O(h3k ).
2
Como queremos hacer hk → 0, podemos asumir que hk < 2/C para la constante que
aparece en la desigualdad anterior, y esto implica que k``k+1 k = O(h3k ).
donde !
s
X
Ki := f tk + ci hk , yk + hk aij Kj , (8.27)
j=1
199
Métodos Numéricos Facultad de Ingenierı́a – IMERL
yk+1 = yk + hk (b1 K1 + b2 K2 ),
donde
K1 = f (tk + c1 hk , yk ) , K2 = f (tk + c2 hk , yk + hk a21 K1 ) .
Por lo tanto, debemos hallar valores de los cinco parámetros b1 , b2 , c1 , c2 , a21 que logren
que el error de truncamiento sea ` k+1 = O(h3k ). Comenzamos, como ya en este punto el
lector podrá sospechar, por el desarrollo de Taylor (8.24):
h2k 00
u(tk+1 ) = yk + hk f (tk , yk ) + u (tk ) + O(h3k ).
2
Para alivianar la notación, pondremos f := f (tk , yk ), ft = ∂f (t , yk ), fy =
∂t k
∂f
(t , yk );
∂y k
las
dos primeras cantidades son vectoriales y la última es matricial. Como
d
u00 (tk ) = f (t, y(t)) t=tk
= ft + fy y0 (tk ) = ft + fy f ,
dt
5
El primer punto es elemental de verificar. La verificación de los dos últimos no es tan directa como
parece.
200
Métodos Numéricos Facultad de Ingenierı́a – IMERL
podemos escribir
h2k h2
u(tk+1 ) = yk + hk f + ft + k fy f + O(h3k ). (8.28)
2 2
Dejamos esta fórmula en pausa y buscamos una expresión análoga para el paso del método
de Runge-Kutta. Para simplificar el análisis, asumimos que c1 = 0, lo que implica K1 = f .
Además, un desarrollo de Taylor para f nos da
K2 = f (tk + c2 hk , yk + hk a21 K1 ) = f + c2 hk ft + fy hk a21 K1 + O(h2k ).
Con estas dos expresiones, obtenemos
yk+1 = yk + hk b1 K1 + hk b2 K2
= yk + hk b1 f + hk b2 [f + c2 hk ft + fy hk a21 f + O(h2k )] (8.29)
= yk + hk (b1 + b2 )f + h2k b2 c2 ft + h2k a21 b2 fy f + O(h3k ).
Para lograr un método cuyo error local de truncamiento sea O(h3k ), basta con comparar los
desarrollos (8.28) y (8.29) e imponer igualdades término a término en los lados derechos:
b1 + b2 = 1
b2 c2 = 12 . (8.30)
1
a21 b2 = 2
Enfatizamos que cualquier elección de b1 , b2 , c2 , a21 que sea solución de este sistema da
lugar a un método de Runge-Kutta explı́cito de segundo orden. Es sencillo observar que el
sistema de arriba tiene infinitas soluciones, que se pueden parametrizar como b2 = α 6= 0,
1
b1 = 1 − α, c2 = a21 = 2α . Mencionamos a continuación dos ejemplos clásicos.
201
Métodos Numéricos Facultad de Ingenierı́a – IMERL
(hλ)2
Por lo tanto, necesitamos que valga |1 + hλ + 2
| < 1. Operando elementalmente,
concluimos la condición de estabilidad absoluta
2
h< .
|λ|
Es interesante observar que esta condición es la misma que la que obtuvimos para el
método de Euler (8.22). Remarcamos, sin embargo, que si el análisis se hace sobre el
plano complejo, en el espı́ritu de la Observación 8.3.4, entonces se obtiene una región de
estabilidad que contiene estrictamente a la del método de Euler.
K1 = f (tk , yk ) ,
hk hk
K2 = f tk + , yk + K1 ,
2 2
yk+1 = yk + hk K2 ,
o, en forma compacta,
hk hk
yk+1 = yk + hk f tk + , yk + f (tk , yk ) .
2 2
2
h< .
|λ|
202
Métodos Numéricos Facultad de Ingenierı́a – IMERL
en tomar
K1 = f (tk , yk ),
hk hk
K2 = f tk + , yk + K1 ,
2 2
hk hk
K3 = f tk + , yk + K2 , (8.33)
2 2
K4 = f (tk + hk , yk + hk K3 ) ,
hk
yk+1 = yk + (K1 + 2K2 + 2K3 + K4 ) .
6
Este es un método de cuarto orden y es condicionalmente absolutamente estable.
Figura 8.5: Solución a (8.11) en el intervalo [0, 5], para a = 1, y0 = 0,1, con paso constante
h = 1. La solución exacta se representa con la curva en negro, mientras que los las
aproximaciones por los métodos de Euler hacia adelante (azul), de Heun (rojo) y de
Runge-Kutta de 4 pasos (verde) se representan con puntos.
0,1e5
y(5) = ≈ 0,9428.
0,1e5 + 1 − 0,1
203
Métodos Numéricos Facultad de Ingenierı́a – IMERL
La siguiente tabla muestra la evolución de los errores en T = 5 para los tres métodos al
tomar pasos de longitud constante e ir reduciendo la longitud del paso a la mitad.
Es interesante observar cómo evolucionan los errores: a grosso modo, cuando se divide h
por la mitad, para el método de Euler el error se reduce a la mitad, para el de Heun en un
factor de 4, y para el de Runge-Kutta de 4 etapas en un factor de 16. Esto es coherente
con los hechos de que el método de Euler es de primer orden, el de Heun de es de segundo
orden, y el de Runge-Kutta de 4 etapas es de cuarto orden.
En general, usar métodos de orden más alto tiene sentido en la medida en que la solución
buscada sea suave y sus derivadas de orden alto no crezcan demasiado rápidamente. 4
Ejemplo 8.6.1 (otra forma de entender el método de Heun). Dado el problema (8.6),
hacemos un par predictor-corrector usando el método de Euler hacia adelante (8.9) como
204
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Enfatizamos que la diferencia entre la segunda fórmula arriba y (8.23) está en que aquı́,
(P )
del lado derecho no aparece el valor desconocido yk+1 sino el valor ya computado yk+1 .
Hemos obtenido el iterado
f (tk , yk ) + f (tk+1 , yk + hk f (tk , yk ))
yk+1 := yk + hk ,
2
que es exactamente igual a (8.31): el método de Heun es el esquema predictor-corrector
que surge al combinar los métodos de Euler hacia adelante con el del trapecio. 4
Ejemplo 8.6.2. Construimos el par predictor-corrector usando el método de Euler hacia
adelante (8.9) como predictor y el de Euler hacia atrás (8.10) como corrector. Para cada
k ≥ 0, dado yk , procedemos en dos etapas:
(P )
yk+1 := yk + hk f (tk , yk ) (paso predictor),
(P )
yk+1 := yk + hk f (tk+1 , yk+1 ) (paso corrector).
(P )
Podemos reemplazar el valor de yk+1 en la segunda fórmula, lo que da lugar a
No nos hemos encontrado con este método directamente, pero podemos interpretarlo
como un método de Runge-Kutta de dos etapas, en las que en (8.26)–(8.27) se toma
s = 2, c1 = 0, c2 = 1, a11 = a12 = a22 = 0, a21 = 1, b1 = 0, b2 = 1. Como este conjunto de
parámetros no satisface (8.30), este método no es de segundo orden. Dejamos al lector la
tarea de verificar que es de primer orden. Este método no tiene mayor interés práctico. 4
8.7. Solvers
Tanto Matlab como Octave tienen incorporados varios solvers, que permiten obtener su-
cesiones {(tk , yk )} utilizando distintos métodos. Naturalmente, distintos métodos pueden
tener diferentes ventajas o desventajas, ser más o menos precisos, o más o menos deman-
dantes computacionalmente. La elección de qué método usar puede depender de nuestros
requerimientos como también de algunas caracterı́sticas del problema que queremos re-
solver.
Los solvers de ecuaciones diferenciales ordinarias son funciones de Octave/Matlab cuyo
nombre es del tipo odeNNX. La sigla ode proviene del inglés ordinary differential equation,
205
Métodos Numéricos Facultad de Ingenierı́a – IMERL
NN es un par de números que indican el orden del método que el solver utiliza, y X es un
conjunto (posiblemente vacı́o) de caracteres que indica alguna caracterı́stica especial de
dicho método. Todavı́a no hemos definido qué es el orden de un método (ver la Definición
8.3.4 más abajo) pero, por el momento, basta con considerar que el orden es un indicador
de precisión.
Una caracterı́stica importante de estos solvers es que utilizan pasos variables, esto es,
dada una cierta tolerancia de error, intentan tomar los pasos lo mayor posible de modo
que la solución computada sea aceptable para dicha tolerancia y que pueda ser computada
de forma económica. La determinación de los pasos se hace usando estimadores de error:
explicamos algunas ideas sobre el tema en la Sección 8.7.1.
A continuación, damos una lista de algunos solvers que podemos utilizar y una breve
descripción como referencia.
ode23 también usa un par de métodos de Runge-Kutta explı́citos, pero uno de orden
2 y otro de orden 3. En forma similar a ode45, los iterados se computan usando el
método de orden 2 y el de orden 3 se usa principalmente para estimar el error y
determinar el largo del paso. Puede ser más eficiente que ode45 para tolerancias de
error grandes o en problemas moderadamente rı́gidos.
ode15s también es un solver para problemas rı́gidos, pero de mayor orden que
ode23s.
En la lista de arriba aparecieron varios términos que aún no hemos definido. En las si-
guientes secciones vamos a profundizar sobre éstas y otras nociones relacionadas a la
aproximación numérica de ecuaciones diferenciales ordinarias. El siguiente ejemplo mues-
tra cómo se utilizan los solvers en la práctica.
206
Métodos Numéricos Facultad de Ingenierı́a – IMERL
solamente pueden comer conejos. Modelamos este ecosistema con un par de ecuaciones
diferenciales no lineales acopladas,
0
c (t) = 2c(t) − αc(t)z(t),
z 0 (t) = −z(t) + αc(t)z(t),
(8.34)
c(t0 ) = c0 ,
z(t0 ) = z0 .
Aquı́, la variable t representa el tiempo, c(t) (resp. z(t)) es la población de conejos (resp.
zorros) en el instante t, y α > 0 es una constante que indica la probabilidad de encuen-
tro entre zorros y conejos. Se puede probar que las soluciones del problema (8.34) son
periódicas, con un perı́odo que depende de las condiciones iniciales. Esto es, dados c0 , z0 ,
existe un tiempo T tal que
Para resolver numéricamente el sistema (8.34), introducimos la función y(t) = [c(t); z(t)],
y escribimos las ecuaciones compactamente,
2y1 (t) − αy1 (t)y2 (t) c
0
y (t) = f (t, y(t)) := , y(t0 ) = y0 := 0 .
−y2 (t) + αy1 (t)y2 (t) z0
alpha = 0.02;
y0 = [100 10];
f = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
Para resolver el problema computacionalmente usando uno de los solvers que mencionamos
más arriba, basta con llamarlo. Por ejemplo, para usar ode23 para resolver el sistema en
el intervalo de tiempo [0, 10], escribimos
207
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 8.6: Solución al modelo de Lotka-Volterra usando los solvers ode23 (izquierda),
ode45 (centro), ode23s (derecha). En la fila de arriba, graficamos las poblaciones de cone-
jos (azul) y zorros (rojo) en función del tiempo. La fila de abajo representa conjuntamente
las poblaciones de conejos y zorros.
Eventos
En el Ejemplo 8.7.1, las soluciones son periódicas. Si queremos determinar el perı́odo
de estas soluciones, esto es, el menor T que cumple (8.35), podrı́amos intentar hacerlo
visualmente a partir de la Figura 8.6: para los valores de los parámetros que tomamos,
parece que T es un poco mayor que 6. Naturalmente, podrı́amos determinar T de forma
un poco más precisa si usamos los valores de poblaciones obtenidos.
Una caracterı́stica interesante de los solvers son los eventos. Uno puede indicar a la
función que tome determinadas acciones si ocurren situaciones de interés. Como primer
ejemplo, pongamos que buscamos el primer tiempo tal que la población de conejos es
igual a 200. Para ello, creamos una nueva función, que en este caso llamaremos parar:
La primera salida, val, es el valor que queremos que sea cero. Fijando la segunda salida
isterm a 1 se indica que el solver debe terminar una vez que val es igual a cero. Poner
que la tercera salida, dir, sea la matriz nula, indica que el cero que estamos buscando
puede ser alcanzado desde cualquier dirección, desde valores positivos o negativos.
Con esto, modificamos nuestro código para llamar al solver como se describe a continua-
ción.
alpha = 0.02;
y0 = [100 10];
F = @(t,y) [2*y(1) - alpha*y(1)*y(2); -y(2)+alpha*y(1)*y(2)];
opts = odeset(’events’,@parar);
208
Métodos Numéricos Facultad de Ingenierı́a – IMERL
209
Métodos Numéricos Facultad de Ingenierı́a – IMERL
iterado computado y volver a intentar con un paso más corto. En la práctica, sin embargo,
uno no tiene a mano la solución exacta del problema, y se suele reemplazar el error por
una cantidad computable: se usa un paso con un método más preciso (tı́picamente de
mayor orden) como solución aproximada y se usa la discrepancia entre las dos soluciones
computadas como estimador del error.
Los métodos de un paso permiten adaptar las longitudes de los pasos de forma relativa-
mente sencilla, porque uno simplemente puede ir cambiando los hk de un paso al siguiente.
Debido a que la de Runge-Kutta es una familia de métodos de un paso que pueden ser
de orden arbitrario, estos métodos son la base de una buena parte de los solvers de
Matlab/Octave. Aquı́ vamos a comentar algunas ideas clave relacionadas al algoritmo
BS237 , que es la base del solver ode23. Para simplificar la notación, vamos a considerar
un problema de valores iniciales (8.6) escalar, esto es, con n = 1. En el caso vectorial, se
trabaja coordenada a coordenada y se debe reemplazar los valores absolutos por normas.
El algoritmo BS23 usa dos métodos de Runge-Kutta combinados, uno de orden 3 y otro
de orden 2. El método de orden 3 es de 3 etapas, mientras que el de orden 2 utiliza las
mismas 3 etapas y le agrega una cuarta; sin embargo, el algoritmo usa el llamado enfoque
FSAL8 , que implica que en la práctica se tengan que hacer 3 evaluaciones de f en cada
paso.
Concretamente, dadas tolerancias de error relativo rtol y absoluto atol, el iterado yk , y
un paso tentativo hk , para determinar el iterado siguiente yk+1 el algoritmo BS23 comienza
computando
K1 := f (tk , yk ) ,
hk hk
K2 := f tk + , yk + K1 ,
2 2
3hk 3hk
(8.36)
K3 := f tk + , yk + K2 ,
4 4
hk
yk+1 := yk + (2K1 + 3K2 + 4K3 ) .
9
El valor de yk+1 obtenido corresponde a un método de Runge-Kutta de tercer orden. A
continuación, se reutilizan los valores K1 , K2 , K3 ya computados para generar un método
de segundo orden. Para ello, se agrega el valor
K4 := f (tk + hk , yk+1 ) ;
notar que en caso de que el paso que se propone sea aceptado, este valor de K4 es el valor
de K1 del paso siguiente. Se calcula un paso estimado
hk
ỹk+1 := yk + (7K1 + 6K2 + 8K3 + 3K4 ) .
24
7
Bogacki, P. y Shampine, L. F. , A 3(2) pair of Runge–Kutta formulas, Applied Mathematics Letters,
2 (1989), pp. 1–9.
8
First Same As Last, “primero igual al último”.
210
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Los factores de 5 y de 4/5 arriba son para prevenir cambios excesivos de longitud entre
pasos consecutivos, y la potencia 1/3 está relecionada con que el método de Runge-Kutta
utilizado es de orden 3.
8.8. Rigidez
Esta sección trata sobre un concepto sutil e importante en la resolución numérica de
ecuaciones diferenciales ordinarias: la llamada rigidez. Consideremos la siguiente variante
sencilla del problema test (8.21):
(
y 0 (t) = λ(y(t) − g(t)) + g 0 (t),
(8.37)
y(0) = y0 ,
donde g es una función conocida y regular, y λ < 0. Se puede verificar que la solución de
este problema es
y(t) = (y0 − g(0))eλt + g(t).
Notemos que, si λ 0, entonces el primer término del lado derecho arriba es despreciable
cuando t es grande: la solución se acerca a g(t) exponencialmente rápido respecto a t. Sin
embargo, si usamos varios de los métodos que hemos discutido hasta ahora para resolver
este problema, nos vamos a encontrar con algunas dificultades. Para que los métodos
explı́citos sean estables se requiere que los pasos sean extremadamente pequeños.
211
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Para fijar ideas, tomamos λ = −10, g(t) = t, y0 = 1, y resolvemos (8.37) sobre el intervalo
[0, 10]. Si usamos el método de Euler hacia adelante o el de Heun con paso constante,
estamos en serios problemas: para que sean estables vamos a necesitar que los pasos
sean menores9 que h = 2/10. La Figura 8.8 muestra soluciones con paso constante y
usando 51 (estable), 50 (lı́mite de estabilidad) y 49 (inestable) pasos, respectivamente. El
comportamiento que observamos es similar al de la primera fila en la Figura 8.4. A pesar
de que la solución es prácticamente lineal para t & 1, con los métodos de Euler y de Heun
no podemos dar pasos muy largos sin perder el control de las soluciones numéricas.
Figura 8.8: Solución a (8.37) en el intervalo [0, 10], para λ = −10, g(t) = t, y0 = 1, con
paso constante. La primera columna muestra las soluciones con el método de Euler y la
segunda columna con el método de Heun, tomando h = 10/51 (arriba), h = 1/5 (centro),
y h = 10/49 (abajo). La solución exacta se representa con las curvas en rojo, mientras
que las aproximaciones numéricas se muestran en azul.
Para este problema, incluso los solvers ode23 y ode45, que usan métodos de Runge-
Kutta explı́citos, están en dificultades. Si bien la solución de nuestro problema es suave,
no podemos usar pasos de tiempo largos. Para mayor dramatismo, cambiamos el intervalo
de resolución a [0, 100] y tomamos λ = −100 en (8.37). La Figura 8.9 muestra la solución
computada con el solver ode45 con los mismos parámetros pero sobre el intervalo [0, 100],
y un zoom en el intervalo [95, 100]. Para este problema, ode45 tomó 12089 pasos para
9
No vamos a entrar en el detalle de por qué se tiene esta condición de estabilidad, pero mencionamos
que se deduce de un razonamiento como el de la Observación 8.3.3. A los efectos de la discusión sobre
rigidez, basta con tener en cuenta que estos métodos son condicionalmente estables y que por lo tanto
hay una longitud de paso máximo admisible.
212
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 8.9: Solución a (8.37) en el intervalo [0, 100], para λ = −100, g(t) = t, y0 = 1, con
el solver ode45. En la izquierda, vemos la solución en todo el intervalo (se ve como una
curva gruesa debido a la enorme cantidad de puntos que se representan), y a la derecha
el detalle para t ∈ [95, 100].
Estamos ante un problema que parece ser sencillo pero que resulta costoso de resolver con
métodos que no sean incondicionalmente estables. En particular, en virtud de la Obser-
vación 8.3.5, cualquier método explı́cito de paso simple va a darnos dolores de cabeza en
esta situación. A problemas de este tipo se los suele llamar rı́gidos. No hay una definición
rigurosa de rigidez, aunque tı́picamente refiere al hecho de que se busca una solución que
varı́a lentamente (es suave), pero que está rodeada por soluciones que varı́an rápidamente.
Esto hace que, si se le aplican métodos no absolutamente estables con control de pasos, los
pasos deban ser extremadamente cortos para obtener resultados satisfactorios. La rigidez
es una cuestión de eficiencia.
Para tratar con problemas rı́gidos, se suelen usar métodos implı́citos. La letra s en los
nombres de solvers como ode23s u ode15s hace referencia precisamente a que usan al-
goritmos adecuados para problemas rı́gidos. Para comparación, la Figura 8.10 muestra
la solución obtenida con el solver ode23s con los mismos parámetros que usamos para
producir la Figura 8.9. En este caso, el solver solamente toma 47 pasos en todo el intervalo
[0, 100], y de esos 47 pasos, 36 se encuentran concentrados en el intervalo [0, 1].
10
Por defecto, en ode45 la tolerancia de error relativa es 10−3 y la absoluta es 10−6 .
213
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Figura 8.10: Solución a (8.37) en el intervalo [0, 100], para λ = −100, g(t) = t, y0 = 1,
con el solver ode23s.
214
Apéndice A
Preliminares y repaso
rk (h)
Con esta notación, el lı́mite (A.1) se puede reescribir como lı́mh→0 hk
= 0. 4
Definición A.1.1 (polinomio de Taylor). El polinomio involucrado en el Teorema A.1.1
se denomina polinomio de Taylor de orden k de la función f en el punto a y lo
denotamos por Tk f 1 ,
f 00 (a) f (k) (a)
Tk f (x) := f (a) + f 0 (a)(x − a) + (x − a)2 + . . . + (x − a)k .
2! k!
4
Observación A.1.2. Una noción intuitiva sobre el polinomio de Taylor Tk f (x; a) es que es
el polinomio de grado k que “mejor aproxima” a f alrededor del punto a. Concretamente,
usando el Teorema A.1.1 y en particular la fórmula (A.1), tenemos
|f (x) − Tk f (x; a)| = |rk (x − a)| = o((x − a)k ).
216
Métodos Numéricos Facultad de Ingenierı́a – IMERL
El Teorema A.1.1 nos da el polinomio que “mejor aproxima” a una cierta función f
alrededor de un cierto punto a en el sentido de que se cumple (A.1). Sin embargo, ese
teorema no nos indica ninguna forma de expresar el resto rk . Concluimos esta sección con
un resultado útil al respecto3 .
Teorema A.1.2 (forma de Lagrange del resto). Si f ∈ C k+1 , entonces, para todo x ∈ R
existe un punto θx entre a y x (esto es, θx ∈ [a, x] si x > a o θx ∈ [x, a] si x < a) tal que
el resto de Taylor se puede expresar como
f (k+1) (θx )
rk (x − a) = (x − a)k+1 .
(k + 1)!
Observación A.1.3. El Teorema A.1.2 nos asegura que el resto de Taylor se puede escribir
de forma análoga que los demás términos en el desarrollo de Taylor, pero con la derivada
(k + 1)-ésima evaluada en un cierto punto θx entre a y x. Es importante notar que θx
depende tanto de x como de a (y obviamente de la función f ), y que el teorema anterior
no nos brinda su valor exacto. 4
Retomamos brevemente el tema de desarrollos de Taylor en una variable en el Capı́tulo 1
de las notas del curso.
217
Métodos Numéricos Facultad de Ingenierı́a – IMERL
El siguiente es una generalización del Teorema A.1.1 (con k = 1) al caso multivariado: dice
que una función diferenciable se puede aproximar, en cada punto de su dominio, con una
transformación lineal afı́n –su polinomio de Taylor de primer grado en varias variables– y
que esta aproximación tiende a cero (en cada coordenada) más rápido que kx − ak.
kr(x − a)k
f (x) = f (a) + Jf (a)(x − a) + r(x − a), y lı́m = 0.
x→a kx − ak
Observación A.2.1. La condición f (a)f (b) < 0, equivale a decir que f toma valores con
signo distinto en los puntos a y b. Esto es, o bien f (a) > 0 y f (b) < 0, o bien f (a) < 0 y
f (b) > 0. 4
Por otra parte, el Teorema de Rolle es un resultado de valor medio para derivadas: nos
da una condición para que la derivada de una función tenga una raı́z.
218
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Darboux: Para todo λ entre f (a) y f (b), existe c ∈ (a, b) tal que f (c) = λ.
f (b)−f (a)
Lagrange: Si f es derivable en (a, b), existe c ∈ (a, b) tal que f 0 (c) = b−a
.
Av = λv.
Aquı́, K es el cuerpo con el que estamos trabajando. En este curso utilizamos principal-
mente el cuerpo de los reales K = R, aunque por momentos consideraremos el cuerpo de
los complejos K = C. 4
Una forma de hallar los valores propios de una matriz es encontrando las raı́ces de su
polinomio caracterı́stico. Este polinomio se define de la siguiente forma.
Teorema A.3.1. Los valores propios de A son las raı́ces del polinomio caracterı́stico de
A en el cuerpo K. Esto es, λ es un valor propio si y sólo si se cumple pA (λ) = 0.
Una vez que se tienen los valores propios de la matriz A, podemos buscar los vectores
propios asociados a cada valor propio.
219
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Teorema A.3.2. Los vectores propios de A asociados al valor propio λ son todos aquellos
vectores v 6= 0 que son solución del sistema lineal (compatible indeterminado)
Av = λv (⇔ (A − λI) v = 0).
En otras palabras, los vectores propios de A asociados al valor propio λ son aquellos
vectores no nulos que pertenecen al núcleo de la matriz A − λI.
Sλ := {v ∈ Rn : Av = λv}.
A = P DP −1 .
Teorema A.3.3. Una matriz A ∈ Mn (R) es diagonalizable si y sólo si existe una base
de Rn formada por vectores propios de A. En ese caso, la matriz diagonal D tiene a los
valores propios de A en su diagonal, y la matriz P tiene como columnas a la base de
vectores propios.
A.4. Ortogonalidad
Lectura recomendada: Libro rojo de GAL 2, Capı́tulo 4 (“Producto interno y
norma”).
Además de permitirnos definir normas, los productos internos nos permiten hablar de
ángulos entre vectores. En particular, podemos generalizar la noción de perpendicularidad.
220
Métodos Numéricos Facultad de Ingenierı́a – IMERL
4
n
Observación A.4.2. Sea S ⊂ R un subespacio vectorial. Si {w1 , . . . , wk } es una base
ortogonal (pero no ortonormal) de S, entonces podemos normalizarla tomando ui :=
wi
kwi k
= √ wi y tendremos, para v ∈ Rn ,
hwi ,wi i
k k
* + k
X X wi w X hv, wi i
PS (v) = hv, ui iui = v, p p i = wi .
i=1 i=1 hwi , wi i hwi , wi i i=1
hwi , wi i
4
221
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Un resultado clave para nosotros va a ser que todo subespacio S ⊂ Rn tiene una base
ortonormal.
hv2 , w1 i
w2 = v2 − w1 .
hw1 , w1 i
k−1
X hvk , wi i
wk := vk − wi ,
i=1
hwi , wi i
222
Métodos Numéricos Facultad de Ingenierı́a – IMERL
La incógnita es una función y(t) definida y derivable hasta orden k para todo t ∈
I ⊂ R. Nuestro objetivo es hallar esta función.
En este material solamente mencionamos formas para resolver analı́ticamente ciertas ecua-
ciones de primer y segundo orden.
y 0 (t)
Z Z
dt = B(t) dt + C.
A(y(t))
Notar que la igualdad vale a menos de una constante aditiva C ∈ R. Haciendo el cambio
de variable t 7→ y(t) en la integral de la izquierda, tenemos
Z Z
1
dy = B(t) dt + C. (A.3)
A(y)
223
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Esto sugiere el siguiente procedimiento para tratar con (A.4) o (A.5): primero se busca
la solución general del problema homogéneo, y luego se busca una solución particular de
la ecuación no homogénea. Una forma de lograr esto último es mediante el método de
variación de constantes o el método de selección.
224
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Para buscar una solución particular, el método de variación de constantes propone buscar
una solución del tipo
yP (t) = C(t)e−at ,
donde ahora C(t) es una función desconocida. Derivando e imponiendo que yP satisfaga
(A.4), nos encontramos con la condición
C 0 (t)e−at = f (t).
Basta con hallar una función C que cumpla esta ecuación para tener nuestra solución
particular.
C 0 (t)e−2t = et .
Por lo tanto, tiene que ser C 0 (t) = e3t e integrando encontramos la solución particular
3t t
yP (t) = e3 e−2t = e3 . Concluimos que todas las soluciones de (A.6) se pueden escribir
como
et
y(t) = Ce−2t + , con C ∈ R.
3
4
Para hallar la solución general del problema homogéneo asociado a la ecuación de se-
gundo orden (A.5), se considera la ecuación caracterı́stica
p(λ) := λ2 + aλ + b. (A.7)
Conociendo las raı́ces de este polinomio asociado, se puede hallar la solución general de
la ecuación homogénea.
225
Métodos Numéricos Facultad de Ingenierı́a – IMERL
Teorema A.5.2 (Solución general homogénea). Sean λ1 y λ2 las raı́ces del polinomio
(A.7). Entonces, la solución general de la ecuación homogénea asociada a (A.5) es
C 1 e λ1 t + C 2 e λ2 t , λ1 6= λ2 , λ1 , λ2 ∈ R,
yH (t) = (C1 + C2 t) eλ1 t , λ1 = λ2 ∈ R,
αt
e (C1 cos (βt) + C2 sin (βt)) , λ1 = α + iβ = λ2 ∈ C.
226
Bibliografı́a
[GVL13] G. Golub and C. Van Loan. Matrix computations. Johns Hopkins University
Press, 2013.
[Mol04] C. Moler. Numerical computing with MATLAB. SIAM, 2004. Capı́tulos dispo-
nibles gratuitamente en [Link]
[Mol20] C. Moler. How the SVD saves the universe? Video de charla. Dis-
ponible en [Link]
Matematick%C3%A9probl%C3%A9mynematematik%C5%AF, 2020.
[Ngu21] Trung Nguyen. The convergence of the regula falsi method. arXiv preprint
arXiv:2109.03523, 2021.
[QS06] A. Quarteroni and F. Saleri. Cálculo cientı́fico con MATLAB y Octave. Springer,
2006.
227