ANÁLISIS

NUMERICO I
 Análisis Numérico / Análisis Numérico I

Clase 1
En esta clase veremos algunas definiciones básicas, presentaremos algunos problemas que

21
estudiaremos durante el curso y repasaremos algunos conceptos matemáticos que serán ne-
cesarios para las clases siguientes.

20
Preliminares y definiciones básicas

Los modelos matemáticos son la herramienta básica para resolver problemas científicos.

F
Generalmente, tales problemas se originan en situaciones o “problemas de la vida real” de

A
las más variadas áreas o disciplinas. Características propias del problema a estudiar permiten
definir ecuaciones y/o inecuaciones que ayudan a formular el modelo matemático.

“vida real”
M
problema de la modelo
matemático
A
-F

Es fácil imaginar que, frecuentemente, en problemas complejos se realizan algunas simplifi-

caciones al construir el modelo. De allí que el modelo matemático no es necesariamente una
descripción exacta de la realidad y que las respuestas que se obtienen deben ser chequeadas
y comparadas con resultados experimentales.
o

Así, un modelo matemático es una buena descripción del problema a estudiar, pero en ge-
neral no siempre da una respuesta directa al problema considerado. Una primera posibilidad
ic

podría ser realizar muchas más simplificaciones en el modelo de manera de poder obtener
una solución analítica. La desventaja de esto es que el modelo matemático podría diferir
ér

mucho del problema real. En lugar de esto, se pueden realizar aproximaciones numéricas
del problema, que si bien pueden introducir errores, es posible estudiar y estimar cómo estas
um

aproximaciones afectan a la precisión de la solución.

modelo matemático 1 modelo matemático 2

.N

problema matemático problema matemático

A

aproximaciones

solución analítica problema numérico

solución numérica

Definición 1 Un problema numérico es una clara y nada ambigua descripción de la cone-

xión funcional entre datos de entrada (variables independientes del problema o input) y los
datos de salida (resultados deseados o output). Estos datos, input y output, consisten en un
conjunto finito de cantidades reales.

1
 Para resolver un problema numérico se disponen, además de la teoría matemática, de 2 he-
rramientas fundamentales: las computadoras y los algoritmos.

Definición 2 Un algoritmo para un problema numérico es una completa descripción de un

número finito de pasos con operaciones bien definidas, y sin ambigüedades, a través de las
cuales una lista de datos de entrada se convierte en una lista de datos de salida.

21
El objetivo del Análisis Numérico es formular y estudiar métodos numéricos y algoritmos
para obtener la solución de problemas provenientes de la vida real y/o de diversas áreas,

20
modelizados matemáticamente.

F
Temario y algunos problemas:

A
1. Teoría de errores. Fuentes de error, aritmética finita, artimética del computador.

M
A
-F
o

2. Ecuaciones no lineales.
R R
Dada f : → hallar x∗ solución de f (x) = 0.
ic
ér
um
.N

3. Interpolación polinomial.
A

2
 4. Teoría de mejor aproximación. Método de cuadrados mínimos.

21
20
5. Integración numérica.

F
A
M
A
-F

6. Resolución numérica de sistemas lineales. Resolver numéricamente sistemas linea-

les de la forma Ax = b.
o

7. Programación lineal.
ic
ér
um
.N
A

Preliminares matemáticos

Teorema 1 (Valor intermedio para funciones continuas). Sea f continua en [a, b]. Sea d
entre f (a) y f (b) entonces existe c ∈ [a, b] tal que f (c) = d.

Teorema 2 (Valor medio). Sea f continua en [a, b] y derivable en (a, b). Entonces para
todo par x, c ∈ [a, b] se cumple que
f (x) − f (c)
= f 0 (ξ ), para algún ξ entre x y c.
x−c
Esto dice que f (x) = f (c) + f 0 (ξ )(x − c).

3
 Teorema 3 (Taylor). Si f ∈ C(n) [a, b] y existe f (n+1) (a, b) entonces para todo par x, c ∈ [a, b]
se tiene que
n
1
f (x) = ∑ f (k) (c)(x − c)k + En (x),
k=0 k!
donde
1
En (x) = f (n+1) (ξ )(x − c)n+1 , para algún ξ entre x y c.
(n + 1)!

21
Observación: tomando y = c, (x − c) = h y por lo tanto x = y + h, entonces

20
h2 00 hn hn+1 (n+1)
f (y + h) = f (y) + h f 0 (y) + f (y) + · · · + f (n) (y) + f (ξ ),
2 n! (n + 1)!
para algún ξ entre y e (y + h).

F
A
Teorema 4 (Taylor con resto integral). Si f ∈ C(n) [a, b] y existe f ()n+1 (a, b) entonces para
todo par x, c ∈ [a, b] se tiene que

M
f (x) =
n
1
∑ k! f (k)(c)(x − c)k + Rn(x),
k=0
A
donde Z x
1
f (n+1) (t)(x − t)n dt.
-F

Rn (x) =
n! c

Órdenes de convergencia
o

Al resolver un problema numérico generalmente no es usual obtener una solución en for-

ic

ma directa o con una fórmula cerrada. Por el contrario, se suele obtener una sucesión de
aproximaciones cuya precisión aumenta progresivamente. De allí que es muy importante el
ér

concepto de convergencia de sucesiones y velocidad de esta convergencia.

Consideremos, por ejemplo, la sucesión dada por xn = 3−n para n ∈ N. Es claro que
um

1
xn = 3−n = −→ 0 cuando n −→ ∞.
3n
Luego
.N

|xn+1 − 0| 3−(n+1) 3n 1
= −n
= n+1
−→ = C.
|xn − 0| 3 3 3
En este caso, cuando C es menor que 1, se dice la convergencia es lineal. En breve formali-
A

zaremos esta definición. Notar que en este ejemplo particular se tiene que
1
|xn+1 − 0| = |xn − 0| para todo n.
3

1
Ahora consideremos la sucesión dada por xn = n! , la que claramente converge a 0. Luego,
|xn+1 − 0| n! 1
= = −→ 0, cuando n → ∞.
|xn − 0| (n + 1)! n + 1
En este caso se dice que la convergencia es superlineal. Veamos formalmente estas defini-
ciones.

4
 Definición 3 Sea {xn } una sucesión de números reales que converge a x∗ .
Se dice que la sucesión {xn } tiene tasa de convergencia (al menos) lineal si existe una
N
constante c tal que 0 < c < 1 y un N ∈ tal que
|xn+1 − x∗ | ≤ c|xn − x∗ | para todo n ≥ N.

Se dice que la tasa de convergencia es (al menos) superlineal si existe una sucesión {εn }
N

21
que converge a 0 y un N ∈ tal que
|xn+1 − x∗ | ≤ εn |xn − x∗ | para todo n ≥ N.

20
Se dice que la tasa de convergencia es (al menos) cuadrática si existe una constante positiva
N
c y un N ∈ tal que
|xn+1 − x∗ | ≤ c|xn − x∗ |2 n ≥ N.

F
para todo

A
Notación O grande y o chica

M
Veremos ahora una notación usual para comparar sucesiones y funciones.
Sean {xn } y {αn } dos sucesiones distintas. Se dice que
A
xn = O(αn )
-F

si existen una constante C > 0 y r ∈ N tal que |xn| ≤ C|αn| para todo n ≥ r.
Se dice que
xn = o(αn )
N tal que |xn| ≤ εn|αn| para
o

si existe una sucesión {εn } que converge a 0, con εn ≥ 0 y un r ∈

todo n ≥ r. Intuitivamente, esto dice que lı́m (xn /αn ) = 0.
ic

n→∞
Ejemplo 1:
ér

 
n+1 1
=O .
n2 n
um

Si
n+1 1 n+1
2
≤C ⇒ ≤ C,
n n n
por lo tanto basta tomar C = 2 y r = 1.
.N

Ejemplo 2:  
1 1
=o .
n ln n n
Si
A

1 1
≤ εn ,
n ln n n
basta tomar εn = 1/ ln n.

Esta notación también se puede usar para comparar funciones. Se dice que
f (x) = O(g(x)) cuando x→∞
si existen una constante C > 0 y r ∈ R tal que | f (x)| ≤ C|g(x)| para todo x ≥ r.
Análogamente, se dice que f (x) = o(g(x)) cuando x → ∞ si lı́m ( f (x)/g(x)) = 0.
x→∞

5
 Ejemplo 3: √ r r
p x 2 +1 x 2 +1 1
x2 + 1 = O(x) si x → ∞, pues = 2
= 1 + 2 ≤ C, luego basta tomar
x x x
C = 2 y r = 1.
Ejemplo 4:
1 1
=o , pues

21
x 2 x
1
x2 x 1
1
= 2
= → 0,
x x

20
x
si x → ∞.
Por último, consideremos el caso de comparación de dos funciones cuando x → x∗ .

F
Se dice que
f (x) = O(g(x)) cuando x → x∗

A
si existen una constante C > 0 y un entorno alrededor de x∗ tal que | f (x)| ≤ C|g(x)| para todo
x en ese entorno.

M
Análogamente, se dice que f (x) = o(g(x)) cuando x → x∗ si lı́m ( f (x)/g(x)) = 0.
x→x∗
A
Ejemplo 5: (Ejercicio)
Ver que
-F

x3
senx = x − + O(x5 ), si x → 0.
6
o
ic
ér
um
.N
A

6
 Clase 2

Algoritmo de Horner: evaluación de polinomios

El algoritmo de Horner es un algoritmo reconocido por su eficiencia para evaluar polinomios

utilizando un número mínimo de operaciones (sumas y productos).

21
Para fijar ideas veamos un ejemplo. Consideremos
p(x) = 2 + 4x − 5x2 + 2x3 − 6x4 + 8x5 + 10x6

20
Primero notemos que, dado x, el método usual para evaluar xk requiere (k − 1) productos:
xk = x| · .{z
. . · }x,

F
k
por ejemplo, si k = 2, se requiere 1 producto, si k = 3, se requieren 2 productos, etc.

A
Método 1: evaluar las potencias de la manera usual, multiplicar por la constante respectiva
y sumar los monomios. Es fácil ver que se requieren:

0+1+2+3+4+5+6 = ∑ i =
6
M 6(6 + 1)
2
= 21 productos, y por lo tanto, se tienen
A
i=1

# sumas = 6
-F

# productos = 21
Método 2: la idea consiste en evaluar las potencias de x en forma sucesiva:
o

x2 = x ∗ x, x3 = x ∗ x2 , x4 = x ∗ x3 , x5 = x ∗ x4 , x6 = x ∗ x5 ,
ic

teniendo en cuenta que cada potencia de x se debe multiplicar por un coeficiente, se tienen
0 + 1 + 2 + 2 + 2 + 2 + 2 = 11 productos, y por lo tanto
ér

# sumas = 6
# productos = 11
um

Método 3: conocido como método de Horner o multiplicación encajada. La idea consiste en

reescribir convenientemente el polinomio p(x) de modo de reducir el número de productos:

p(x) = 2 + x(4 + x(−5 + x(2 + x(−6 + x(8 + x · 10))))).

.N

Es fácil ver que ahora se requieren:

# sumas = 6
A

# productos = 6
Si el grado de p(x) es n, se requieren n(n + 1)/2 productos en el método 1, 2n − 1 en el
método 2 y sólo n en el método 3.
Si p(x) = a0 + a1 x + a2 x2 + · · · + an xn , con an 6= 0, la evaluación de p(x) en x = z se realiza
con los siguientes pasos:
bn−1 = an
bn−2 = an−1 + z ∗ bn−1
..
.
b0 = a1 + z ∗ b1
p(z) = a0 + z ∗ b0 .

1
 En forma más compacta, podemos escribir:

Algoritmo de Horner (multiplicación encajada)

Dados el polinomio p(x), de grado n, con coeficientes ai , para i = 0, . . . n, con an 6= 0 y un

número real z en el que se desea evaluar p(x).

21
input n; ai , i = 0, . . . , n; z
bn−1 ← an

20
for k = n − 1 to 0 step −1, do
bk−1 ← ak + z ∗ bk

F
end do

A
output bi , i = −1, . . . , n − 1
end
Notar que b−1 = p(z) M
A
-F

Fuentes de error y aritmética de la computadora

Una de las tareas más importantes en Análisis Numérico es estimar la precisión del resultado
en un cálculo numérico.
o

Los resultados numéricos están influenciados por diferentes tipos de errores. Algunos de
ic

estos pueden eliminarse, en otros se puede reducir su influencia, y otros son inevitables y
nada se puede hacer.
ér

Principales fuentes de error.

um

1. Errores en los datos de entrada (frecuentemente inevitables). Los datos de entrada

pueden ser el resultado de mediciones experimentales que podrían estar afectadas a
errores sistemáticos debido al equipamiento utilizado. También, se tienen los errores
.N

que se producen al representar un número real irracional con un número finito de

dígitos.
A

2. Errores de redondeo. Aparecen cuando los cálculos se realizan usando un número

finito de dígitos.

3. Errores de truncamiento. Aparecen cuando un proceso infinito es reemplazado por

un proceso finito. Ejemplos: aproximar una serie por una suma parcial, aproximar una
función por un polinomio (Taylor), aproximar una derivada por un cociente incremen-
tal, etc.

4. Errores humanos. Son frecuentes y a veces difíciles de detectar. Ejemplos: errores en

la formulación del problema, en los cálculos “a mano”, al escribir un programa en la
computadora, etc.

2
 Errores absolutos y relativos

Antes de analizar que ocurre al representar números en una computadoras y se realizan ope-
raciones, veremos algunos conceptos básicos que sirven en un contexto más general.

Definición 1 Cuando un número real r (valor exacto) es aproximado por otro número r̄, se
define el error por r − r̄. Llamaremos, respectivamente, a

21
Error absoluto: ∆r = |r − r̄|,
r − r̄ ∆r

20
Error relativo: δr = = .
r |r|
También se llama error (relativo) porcentual al producto 100 ∗ δ r.

F
En general, el error relativo y el error porcentual son más útiles que el error absoluto, porque
da una idea del error cometido relativo a la magnitud de la cantidad que se está considerando.

A
En términos prácticos no se conocen exactamente los valores de los errores absolutos y

M
relativos sino que se tienen cotas de
√ ellos. Siempre se trata que las cotas sean lo más ajustadas
posible. Así, por ejemplo, si r = 2 y r̄ = 1.414, entonces
√
A
∆r = |r − r̄| = | 2 − 1.414| = |0.0002135...| ≤ 0.00022
∆r |0.0002135...|
-F

δr = = √ ≤ 0.00016
|r| 2
y el error porcentual es 0.016 %.
Las siguientes notaciones son equivalentes:
o

r̄ = 1.414, ∆r ≤ 0.22 10−3

ic

r = 1.414 ± 0.22 10−3

ér

1.41378 ≤ r ≤ 1.41422
um

Redondeo y truncado

Existen dos maneras de escribir un número real reduciendo el número de dígitos: redondeo
y truncado.
.N

Redondeo. Veamos algunos ejemplos de redondeo a 3 dígitos (decimales):

0.774432 → 0.774
A

1.23767 → 1.238
0.3225 → 0.323

En general, la aproximación por redondeo a n dígitos decimales r̃ de un número r es un

número de n dígitos decimales (después del punto decimal) que coinciden con los de r si el
dígito (n + 1) de r es 0, 1, 2, 3 o 4. Por otro lado, si el dígito (n + 1) de r es 5, 6, 7, 8 o 9,
entonces para definir r̃ se suma una unidad al n-ésimo dígito de r. Así se cumple que:
1 −n
|r − r̃| ≤ 10 . (1)
2
Para fijar ideas veamos un ejemplo muy sencillo, con n = 1:

3
 si r = 0.11 entonces r̃ = 0.1 y |r − r̃| = 0.01 ≤ 0.05 = 5 10−2 = 12 10−1 ;

si r = 0.16 entonces r̃ = 0.2 y |r − r̃| = 0.04 ≤ 0.05 = 5 10−2 = 12 10−1 .

Truncado. Veamos algunos ejemplos de truncado a 3 dígitos (decimales):

0.774432 → 0.774

21
1.23767 → 1.237
0.3225 → 0.322

20
En general, la aproximación por truncado (o truncamiento) a n dígitos decimales r̂ de un
número r es un número de n dígitos decimales (después del punto decimal) que coinciden
con los n primeros dígitos de r. Así se cumple que:

F
|r − r̂| ≤ 10−n , (2)

A
Para fijar ideas veamos un ejemplo muy sencillo, con n = 1:

M
si r = 0.11 entonces r̂ = 0.1 y |r − r̂| = 0.01 ≤ 0.1 = 10−1 ;
A
si r = 0.19 entonces r̂ = 0.1 y |r − r̂| = 0.09 ≤ 0.1 = 10−1 .
-F

En general, las computadoras trabajan con el sistema de redondeo y no el truncado.

Es importante observar que, si bien el error absoluto que introduce el redondeo depende de
la magnitud del número, el error relativo, que es el más significativo, por ser independiente
o

de la magnitud, está controlado en términos de la cantidad de dígitos con la que se trabaja.

Así vamos a definir el número de dígitos significativos en términos del error relativo.
ic

Definición 2 El número r̄ aproxima a r con m dígitos significativos si

ér

∆r
δr = ≤ 5 10−m .
um

|r|

Esto dice que el error relativo es del orden de 10−m .

.N

Así, en los ejemplos anteriores tenemos

r r̃ ∆r δr díg. signif.
0.774432 0.774 0.00043 0.00056(< 5 10−3 ) 3
A

−4
1.23767 1.238 0.00033 0.00027(< 5 10 ) 4
0.3225 0.322 0.0005 0.0015(< 5 10 ) −3 3

Errores en las operaciones

Analizaremos los errores que se introducen en las operaciones básicas (suma, resta, multi-
plicación y división). Más específicamente, nos centraremos en las dos primeras.
Sean x1 , x2 ∈ R, y x¯1, x¯2 aproximaciones de x1 y x2 respectivamente.
Sean y = x1 + x2 , ȳ = x¯1 + x¯2 .

4
 El error en la operación suma está dado por:

y − ȳ = (x1 + x2 ) − (x¯1 + x¯2 ) = (x1 − x¯1 ) + (x2 − x¯2 ),

por lo tanto el error absoluto

∆y = |y − ȳ| ≤ |x1 − x¯1 | + |x2 − x¯2 |

21
∆y ≤ ∆x1 + ∆x2

y el error relativo
∆y ∆x1 + ∆x2

20
δy = ≤ .
|y| |x1 + x2 |
n n
En general, si y = ∑ xi entonces ∆y ≤ ∑ ∆xi .

F
i=1 i=1

A
El caso de la resta es similar. Sean y = x1 − x2 , ȳ = x¯1 − x¯2 . Entonces el error absoluto se
obtiene haciendo

M
∆y = |y − ȳ| = |(x1 − x2 ) − (x¯1 − x¯2 )| = |(x1 − x¯1 ) − (x2 − x¯2 )|
∆y ≤ ∆x1 + ∆x2
A
y el error relativo
-F

∆y ∆x1 + ∆x2
δy = ≤ .
|y| |x1 − x2 |
o

Es posible estimar cotas de errores para la multiplicación y la división definiendo y = x1 ∗

x2 , ȳ = x¯1 ∗ x¯2 , z = x1 /x2 , z̄ = x¯1 ∗ x¯2 . Se puede deducir que:
ic

∆y ∆x1 ∆x2
∆y / |x2 |∆x1 + |x1 |∆x2 δy = / +
ér

|y| |x1 | |x2 |

y que
um

1 |x1 | ∆z ∆x1 ∆x2

∆z / ∆x1 + 2 ∆x2 δz = / +
|x2 | |x2 | |z| |x1 | |x2 |
.N
A

5
 Clase 3

Cancelación de dígitos significativos

Recordemos la definición de dígitos significativos: el número a aproxima al número real a

con r dígitos significativos si

21
∆a 1
≤ 5 10−r = 101−r .
|a| 2

20
Un efecto no deseable en algoritmos numéricos es la gran cancelación de dígitos significati-
vos que se produce en la resta de números próximos. Para fijar ideas veamos un ejemplo.

F
Sean x1 = 10.123455 ± 0.5 10−6 y x2 = 10.123789 ± 0.5 10−6 .
x1 y x2 tienen error absoluto menor o igual a 0.5 10−6 y error relativo menor a 0.5 10−7 =

A
5 10−8 , esto significa que ambos tienen 8 dígitos significativos.

M
Ahora bien, la resta y = x1 − x2 = −0.000334 ± 10−6 , tiene un error absoluto pequeño, sin
embargo el error relativo
A
∆y 10−6
≤ < 3 10−3 < 5 10−3 ,
|y| 0.000334
-F

por lo tanto la resta y tiene sólo 3 dígitos significativos.

Por lo tanto, es recomendable evitar restas de números próximos, siempre que sea posible.
o
ic

Representación de números en una computadora

ér

El ser humano está acostumbrado a utilizar un sistema de numeración decimal, el cual es un

sistema posicional con base β = 10. La mayoría de las computadoras usa internamente la
base β igual a 2 o 16.
um

Definición 1 sea β ∈ N, β ≥ 2, todo número real r puede ser escrito en la forma:

(±dn dn−1 . . . d2 d1 d0 .d−1 d−2 . . . )β
.N

donde dn , dn−1 , . . . d0 , d−1 . . . son números naturales entre 0 y (β − 1). El valor del número r
es:
±dn β n + dn−1 β n−1 + · · · + d2 β 2 + d1 β 1 + d0 β 0 + d−1 β −1 + d−2 β −2 + . . .
A

Ejemplos:

1. (760)8 = 7 82 + 6 81 + 0 80 = (496)10

2. (101.101)2 = 1 22 + 0 21 + 1 20 + 1 2−1 + 0 2−2 + 1 2−3 = (5.625)10

1
3. (0.333 . . . )10 = 3 10−1 + 3 10−2 + · · · =
3

1
 4. (0.1)10 = (0.0001100110011 . . . )2

Notar que en el último ejemplo (0.1)10 no tiene representación binaria finita!

Observaciones:

1. la mayoría de los números reales no pueden ser representados exactamente en cual-

21
quier base;
2. aparecen errores de representación cuando un número es convertido de un sistema de

20
numeración a otro;
3. aparecen errores debido a que la computadora usa aritmética finita.

F
¿cómo se representan los números en una computadora?

A
Básicamente, existen dos sistemas de representación de números en una computadora:

sistema de punto fijo, M

A
sistema de punto flotante.
-F

Sistema de punto fijo

El primero de ellos es el utilizado por las primeras computadoras (aproximadamente en

1940–1950) y donde los números se representan utilizando una cantidad fija de números
o

enteros y de números fraccionarios. Por ejemplo, si usáramos la base β , (s + 1) dígitos para

ic

la parte entera y t para la parte fraccionaria, tendríamos:

ér

±ds ds−1 . . . d2 d1 d0 . d−1 d−2 . . . d−t ,

donde cada di ∈ {0, . . . , β − 1}. En sistemas contables, aún hoy en día, suele usarse este
um

sistema donde la cantidad de dígitos fraccionarios es t = 2 para representar los centavos.

La principal desventaja de este sistema es que no es posible representar simultáneamente
números reales muy pequeños y muy grandes, sin que la cantidad de dígitos s y t sean
demasiados grandes. Por ejemplo si s = 3 y t = 3, el número más grande y el más que
.N

pequeño que se pueden representar en este sistema son 999.999 y 000.001, respectivamente.
La manera de solucionar este problema es usar la notación científica y esto da origen al otro
sistema.
A

Sistema de punto flotante

Definición 2 Un sistema de punto flotante (β ,t, L,U) es el conjunto de números norma-

lizados en punto flotante en el sistema de numeración con base β , y t dígitos para la parte
fraccionaria, es decir, números de la forma:

x = mβe

donde
m = ±0 . d−1 d−2 . . . d−t

2
 con d−i ∈ {0, . . . , β − 1} para i = 1, . . . ,t, con d−1 6= 0 y L ≤ e ≤ U. Además, β , e y m se
denominan base, exponente y mantisa, respectivamente. Es decir, 1/β ≤ |m| < 1.

Observaciones:

1. aunque el sistema de punto flotante permite representar magnitudes de órdenes muy

21
variados, a diferencia del sistema de punto fijo, también puede ocurrir overflow si
e > U o underflow si e < L;

2. el cero no puede representarse en este sistema de números normalizados.

20
Errores de redondeo en aritmética de punto flotante

F
Al representar números en un sistema de punto flotante (β ,t, L,U) se producen errores de

A
redondeo debido a la precisión limitada. Asumiendo redondeo, estimaremos una cota de los
errores absoluto y relativo.

M
Supongamos que podemos escribir un número real (exacto) en la forma:

x = m β e,
1
A
≤ |m| < 1,
β
-F

donde el exponente e es tal que L ≤ e ≤ U.

Escribimos ahora su representante en el sistema de punto flotante:
1
f l(x) = xr = mr β e , ≤ |m| < 1,
o

β
ic

donde mr es la mantisa que se obtiene redondeando a t dígitos la parte fraccionaria de m.

Entonces, es claro que
ér

1
|mr − m| ≤ β −t ,
2
y por lo tanto, una cota del error del absoluto de representación en x es
um

1
|xr − x| ≤ β −t β e .
2
.N

Para el error relativo, tenemos lo siguiente:

1 −t e
|xr − x| β β 1 −t 1 1−t
≤ 2 e
= β ≤ β ,
|x| |m|β 2|m| 2
A

1
pues si |m| ≥ 1/β entonces ≤ β.
|m|
Luego el error relativo debido al redondeo en la representación en el sistema de punto flotante
está acotado por:
|xr − x| 1 1−t
≤ β = µ,
|x| 2
donde µ se llama unidad de redondeo.
Notar que el error absoluto de representación en punto flotante depende del orden de
la magnitud, en cambio el error relativo no.

3
 ¿Cómo se realiza la suma en aritmética de punto flotante?

Para fijar ideas veremos dos ejemplos con el sistema de punto flotante dado por (β ,t, L,U) =
(10, 4, −9, 9). Sean
x = mx β ex , y = my β ey ,
con x ≥ y. Queremos calcular z = f l(x + y)

21
Ejemplo 1: sean x = 0.1234 100 , y = 0.4567 10−2 , entonces
x + y = 0.1234 100 + 0.4567 10−2

20
= 0.1234 100 + 0.004567 100
= (0.1234 + 0.004567) 100 = 0.12797 100 ,

F
por lo tanto, f l(x + y) = 0.1280 100 .

A
Ejemplo 2: sean x = 0.1234 100 , y = 0.5678 10−5 , entonces

M
x + y = 0.1234 100 + 0.5678 10−5
= 0.1234 100 + 0.000005678 100
A
= (0.1234 + 0.000005678) 100 = 0.123405678 100 ,
-F

por lo tanto, f l(x + y) = 0.1234 100 = x.

Esto ocurre porque el orden de magnitud de x es con respecto a y es muy grande.
o

Observación: algunas propiedades o axiomas de la aritmética infinita dejan de valer en

ic

aritmética de punto flotante. Veamos con un ejemplo que la propiedad asociativa ((a + b) +
c = a + (b + c)) no es válida en un sistema de punto flotante, es decir: f l( f l(a + b) + c) 6=
ér

f l(a + f l(b + c)).

Dado el sistema de punto flotante dado por (β ,t, L,U) = (10, 4, −9, 9) consideremos los
um

números a = 0.9876 104 , b = −0.9880 104 y c = 0.3456 101 . Entonces, por un lado,
f l( f l(a + b) + c) = f l( f l(0.9876 104 − 0.9880 104 ) + c)
= f l( f l(−0.0004 104 ) + c)
.N

= f l(−0.4000 101 + 0.3456 101 )

= f l(−0.0544 101 )
A

= −0.5440 100

Por otro lado,

f l(a + f l(b + c)) = f l(a + f l(−0.9880 104 + 0.0003456 104 ))
= f l(a − f l(0.9876544 104 ))
= f l(0.9876 104 − 0.9877 104 )
= f l(−0.0001 104 )
= −0.1000 101

4
 Observaciones de implementación:

1. dado que en una implementación o programa se realizan muchas operaciones, cada

una con su correspondiente error, es conveniente prestar atención en las operaciones
que se realizan;

2. si x e y son números reales y en una programa se tiene una sentencia del tipo

21
if x == y then . . .

20
es más conveniente reemplazarla por una sentencia del tipo

if (abs(x-y)) < epsilon then . . .

F
para algún valor de epsilon dado por el usuario, puesto que es casi imposible que se

A
verifique la primera sentencia.

M
A
-F
o
ic
ér
um
.N
A

5
 Clase 4 - Solución de ecuaciones no lineales

El problema

Dada f : R → R no lineal, se desea encontrar una solución r de la ecuación

21
f (x) = 0.

R
Idea: comenzando con algún x0 ∈ , generar una sucesión {xk } a través de un algoritmo

20
numérico iterativo, y se espera que tal sucesión converja a r donde f (r) = 0.

Método de bisección

F
Este método se basa fuertemente en el teorema del valor intermedio: si f es continua en [a, b]

A
y si f (a) f (b) < 0, entonces f debe tener una raíz en (a, b).

Idea del método de bisección M

A
a+b
Si f (a) f (b) < 0, se calculan c = y f (c). Sean x0 = c: una aproximación a la raíz r de
2
-F

b−a
f y |e0 | = |x0 − r| ≤ : error de aproximación inicial. Se tienen 3 posibilidades:
2
1. si f (a) f (c) < 0, entonces hay una raíz en el intervalo [a, c]. Reasignamos b ← c y se
Z XÜC
ȈŚ‰ǢžȈržĔȈ,
Ȉ °Þ႒
repite el procedimiento en el nuevo intervalo [a, b].
o

+Ŕಏ
Z XÜC
ȈŚ‰ǢžȈržĔȈ,
Ȉ °Þ႒
ic

П႒
+Ŕಏ
ér

П႒

cBO CMಏ 0,×Û +ŔƴÝಏɵ íಏ +Ýಏ

˜H ,Hƿ ,Gƿ ,ƿ
જБêԢಏ
um

ŏƿ'H, D·ĭ ƿ
cBO CMಏ 0,×Û +ŔƴÝಏɵ íಏ +Ýಏ
˜H ,Hƿ ,Gƿ ,ƿ Figura 1: Caso: f (a) f (c) < 0 જБêԢಏ
ŏƿ'H, D·ĭ ƿ
entonces hay una raíz en el intervalo [c, b]. Reasignamos a ← c y se
2. si f (a) f (c) > 0, +Ŕಏ
repite el procedimiento en el nuevo intervalo [a, b].
.N

+Ŕಏ
۱‫ݩ‬
ƛ‫ݩ‬

+Ý+੣ಏɵ íಏ۱‫ݩ‬
A

ƛ‫ݩ‬
cB O CMಏ 0,ÖÛ
˜H Ĵ,Hƿ , ƿ ,ƿ
DHG,ƿ'H, D· ƿ ƴêಏ
+Ý+੣ಏɵ íಏ
cB O CMಏ 0,ÖÛ
˜H Ĵ,Hƿ , ƿ ,ƿ
DHG,ƿ'H, D· ƿ ƴêಏ

Ɓ€ˆ૵Figura 2: Caso:
ॷ #* ƒ=ॷ/ॷĸɨ _ॷ Ɓ૵ ૵f (a) f (c) > 0 Đ Íù˭
|ॷॷ=/øॷ ॷ / ॷ =Ĝॷ/ॷ ॷ )ॷ/ॷ/ =Lॷĸɨ τॷ  "  ‘
©åࣄ©åɸ૵
r=?)=‘ ?`‘=#E5‘#`‘ E#Z<.‘ #`‘ g=.‘ QZ?<?O#E‘ ?Og.Zp#I ‘  8 ‘©åɸ©åν૵ś Ķ˪ g=.O‘
/ ॷ/ॷ )ॷ ॷॷ=/øॷ R႒ ॷɨ ţ ॷ\/ॷ )=ॷ<ॷ ႒ /ॷ/ =Lॷ/ //ॷ #* ƒņॷRॷ
ॷ Ɓ€ˆ૵
#* ƒ=ॷ/ॷĸɨ
ɏॷ )ॷ Ɓ૵ ૵
&/ॷ:=_ॷ |ॷॷ=/øॷ ॷ / ॷ =Ĝॷ/ॷ
ûॷ/ ॷ )ॷ:=ॷ/ॷ(;ॷ=:
©åࣄ©åɸ૵
ॷ )ॷ/ॷ/
 Ǘॷ Â =ॷ  
#*=Lॷĸɨ
©åɸ©åν૵ś Ķ˪ / ॷ ʻ(Kॷτॷ)ॷ
"‘

3. si f (a) f (c) = 0, entonces f (c)©å૵

r=?)=‘ ?`‘=#E5‘#`‘
)ॷ ॷ = 0Đ˭yÍ˭xĐ˭0 ©åν/૵
E#Z<.‘
<ॷ  h/ॷ = c es Đ Íù˭
#`‘ g=.‘ QZ?<?O#E‘ ?Og.Zp#I ‘  8 ‘€)ॷ‫ܢ‬
la raíz buscada. g=.O‘ ॷ=ॷ )ȇॷ )ॷ
=ॷhÖॷ
/ ॷ/ॷ/ॷ5
( )ॷ ॷॷ=/øॷ R႒ ॷɨ ţ ॷ\/ॷ
) ॷˢ/ ॷ/ॷƒ=ॷ& ॷ/ )=ॷ<ॷ ૵
႒ /ॷ/)ॷ/
ॷॷ )ॷƒ=ॷ/ॷ =Lॷ/ //ॷ
=Lॷĸɨ #* ƒņॷRॷ
_ॷ Ɨ ¦ ॷĪRॷ =ûॷ
Este caso es casi ©å©åӝ૵
ɏॷ
imposible3˪ 
!૵
en la
Ѭpráctica
Rॷ)ॷ &/ॷ:=yå  )/ॷÍ˭
ûॷ/ ॷ )ॷ:=ॷ/ॷ(;ॷ=:
debido  Ǘॷ Â
/ ॷ #* ƒ=ॷ)ॷ &/ॷ
a los errores R =ॷ 
#* / ॷ ʻ(Kॷ
/ Ǘॷ ķL=ûॷ
de redondeo. Por lo )ॷ
 ॷॷI  ॷ
©૵
)ॷ å Ɓ૵
ॷ <ॷ
/Ö'ॷ ˪ Đ˭ Í˭ Đ˭ ©åν/૵
©å૵
 ॷ  ?b‘
h/ॷ
ࣟ 'ॷ /ॷ )ॷ5: =ॷ& €)ॷ‫ܢ‬ =ॷhÖॷ
ॷRॷ= ॷ=ॷǗॷ)ȇॷ
/ ˟ॷ=== ~ ûॷ)ॷॷ
tanto, el criterio de parada no ૵
dependerá
( /ॷ5 ) ॷˢ/ )
`
L˪ 
˪ ©૵åƁ૵ Í <˭
de
ॷ/ॷƒ=ॷ& que f (c)
ॷ/ ॷॷ = 0 ૵
sino de
=ॷ)ॷƒ=ॷ/ॷ que | f (c)| < T OL,
y Û)ॷ/=/Zॷ
=Lॷĸɨ _ॷ Ɨ=/ॷ
¦ ॷĪRॷ =ûॷ

Ѭ
 ॷ (;ॷ) h ॷ
©å©åӝ૵
Rॷ 
!૵
3˪ ɨ Ē åӝί ૵
(;ॷ ûॷyå)ॷॷ
donde T OL es una tolerancia©૵åƁ૵ dada
/Ö'ॷ
ś˪
˪
Í˭ Œ˪
)/ॷ
 ॷ ॷȇ:
ǐpor el usuario. //ॷ ƒ=ॷ)ॷ
ॷ #* ॷ Vf© l &/ॷ
8 3 * ãR
7û˪ / Ǘॷ ķL=ûॷ
J|ॷ| /ॷ   *‫ޓ‬ಏ
?b‘ ࣟ 'ॷ /ॷ )ॷ5: =ॷ& ॷRॷ= / ˟ॷ=== Ǘॷ ~ ûॷ ॷ
ॷॷI R=ॷ
 ॷ #*
=: /ॷRॷ ) ॷ Vf© l  8 3 ã h: =ǗKॷ ~;ॷ( /ॷh  ॷ ॷ ॷ
૵
ŭ
˪
=˪ 1ñ˪ŭʍ1LǮ˪©૵åƁ૵ Í <˭
`
L˪ 
˪ )  ॷ
ɨ Ē åӝί૵ś˪ ǐŒ˪
‫ޔ‬/ॷॷҐ)ॷ
(;ॷ ûॷ )ॷॷ
(;ॷ) =ॷ y Û
7û˪
=/Zॷ =/ॷ h  ॷ
/ॷ ॷVf© l  8 3 * ã J|ॷ| /ॷ   *‫ޓ‬ಏ R=ॷ #*
=: /ॷRॷ ॷȇ: ) ॷ Vf© l  8 3 ã h: =ǗKॷ ~;ॷ( /ॷh  ॷ ॷ ॷ 1
͂ॷ ॷ(Ü /ॷh  ॷ ॷˢ/ ॷ )ॷ= ॷR ॷI  /ॷ” —Ѭ
Ѭ/ Ïॷ  ॷ ॷá˪
‫ޔ‬/ॷॷҐ)ॷŭ
˪
=˪ 1ñ˪ŭʍ1LǮ˪
'H7 =0'— —
 Veamos algunos comentarios de implementación antes de dar el algoritmo.

- Por cuestiones numéricas, en vez de calcular el punto medio haciendo c ← (a + b)/2,

es más conveniente calcular c ← a + (b − a)/2.

- Para determinar el cambio de signo de la función, en vez de analizar si f (a) f (c) < 0,
es más conveniente usar la función sign, y analizar si sign( f (a)) 6= sign( f (b)).

21
- Se utilizan 3 criterios de parada en el algoritmos:

20
1. el número máximo de pasos permitidos (M);
2. el error en la variable es suficientemente pequeño (δ );

F
3. el valor de | f (c)| es suficientemente pequeño (ε).

A
Se pueden construir ejemplos patológicos y simples donde uno de los criterios vale
y el otro no. (Ver Figuras (3 y 4)).Para que el algoritmo sea más robusto se utilizan

M
los tres criterios de parada. De todos modos, siempre es conveniente revisar que se
cumplan las hipótesis para que se pueda aplicar el método y el algoritmo y analizar los
resultados obtenidos.
A
 Z XÜC<
Ȉś ȈƂMĔȈ,
Ȉ °°႒

 Z XÜC<
Ȉś ȈƂMĔȈ,
Ȉ °°႒

-F

ಏ
ಏ
o
ic
ér

ଛ+Ý ৌಏ Ճ‫ ݩڞ ݩ‬/-7 োఋಏ ˅ ë ߈ಏ ୮႒ ྏ႒

ଛ+Ý ৌಏ Ճ‫ ݩڞ ݩ‬/-7 োఋಏ
˅ ë ߈ಏ ୮႒ ྏ႒
cB O CMಏ 
9< ಏ
ಏ
um

cBAH
-ˆN f}y¾ O CMಏB
9<
˄ I1F ±႒
-ˆN f}y¾ AH UDfn
˄ B I1FN¾±႒
UDfn N¾
Figura 3: Ejemplo 1
.N

ࠎ႒
ʡѬ
ࠎ႒
ʡѬ
ಏ
ಏ H႒
<႒H႒
<႒<႒
<႒<႒
<႒H႒
A

<႒H႒
<႒H႒
sѬ
H႒<႒

sѬ
<႒
ú‫ݩ‬
۰‫ݩ‬
ú‫ݩ‬
ë‫ݩ‬ ۰‫ݩ‬ ၲ ȫ ႒ ಏ
ë‫ݩ‬ ၲH႒ ȫ ႒ ಏ
<႒H႒
<႒H႒

ssѬѬ
<႒H႒
<႒
<႒
<႒
ɽcBOCMಏ 
8< <႒
. fS ɽcBOCMಏ
<႒H႒
f}y¾ 1 0
8<
ssѬѬ
I1 > F <႒H႒
. fS f}y¾ 1 0 I1 > F <႒
UDfn‡¾
UDfn‡¾

Figura 4: Ejemplo 2

˪˪vlis
xŽ ×ಏ vlis
yŽ ˪`
˪

i‰
Ѭ
vlis
xŽ ×ಏ vlis
yŽ ˪`
˪

˪i‰

˪
ѬWॷWॷ
˪

˪V˭
V˭
˪

˪ 2
„Ž
„Ž
 „,:Ž
„,:Ž
ŠŠ   „Ž
  „Ž
 Algoritmo de bisección

Dados los siguientes datos de entrada y parámetros algorítmicos: a y b extremos izquierdo y

derecho del intervalo, M el máximo número de pasos (iteraciones) permitidas, δ la tolerancia
para el error e (en la variable x) y ε la tolerancia para los valores funcionales.
input a, b, M, δ , ε

21
u ← f (a)
v ← f (b)

20
e ← b−a
output a, b, u, v

F
if sign(u) = sign(v) then STOP (1)
for k = 1, 2, . . . , M do

A
e ← e/2
c ← a+e
w ← f (c)
M
A
output k, c, w, e
-F

if |e| < δ or |w| < ε then STOP (2)

if sign(w) 6= sign(u) then
b←c
o

v←w
ic

else
ér

a←c
u←w
um

endif
end
.N

Observaciones:
1. En el algoritmo aparecen dos paradas (STOP). La primera detecta que los signos en los
A

extremos del intervalo inicial son iguales y por lo tanto no puede continuar el algoritmo.
Notar que el algoritmo no puede chequear continuidad (ver Figura (4)). La segunda parada
detecta que alguno de los criterios de parada adoptados se cumple.
2. Se podría usar otro criterio de parada basado en que el error relativo sea pequeño en lugar
de considerar el error absoluto.
3. El algoritmo encuentra una raíz de f en [a, b], aunque pueden existir varias raíces. Para
localizar alguna en particular se debe elegir el intervalo inicial cuidadosamente.

3
 El siguiente teorema da un resultado de convergencia del método. Por un lado muestra que
el método es global, en el sentido que si se cumplen las hipótesis del teorema del valor
intermedio, el algoritmo encuentra una solución independiente del tamaño del intervalo. Por
otro lado, de la demostración se puede deducir que la sucesión generada por este algoritmo
converge linealmente.
Teorema 1. Si [a0 , b0 ], [a1 , b1 ], . . . , [an , bn ], . . . denotan los sucesivos intervalos en el método

21
de bisección, entonces existen los límites: lı́m an y lı́m bn , son iguales y representan una
n→∞ n→∞
1 1
raíz de f . Si cn = (an + bn ) y r = lı́m cn , entonces |r − cn | ≤ n+1 (b0 − a0 ).
2 n→∞ 2

20
Demostración. Si [a0 , b0 ], [a1 , b1 ], . . . , [an , bn ], . . . son los intervalos generados en el método
de bisección, se tiene que
a0 ≤ a1 ≤ a2 ≤ · · · ≤ b0

F
b0 ≥ b1 ≥ b2 ≥ · · · ≥ a0 .

A
Luego {an } es una sucesión creciente y acotada superiormente, entonces {an } es convergen-
te. Análogamente, {bn } es una sucesión decreciente y acotada inferiormente, por lo tanto
también es convergente. Además,
M 1
A
bn+1 − an+1 = (bn − an ), n ≥ 0,
2
-F

y aplicando esta ecuación repetidamente, se obtiene que

1
bn − an = (b0 − a0 ), n ≥ 0.
2n
o

Luego
ic

1 1
lı́m bn − lı́m an = lı́m (bn − an ) = lı́m (b0 − a0 ) = (b0 − a0 ) lı́m = 0.
n→∞ n→∞ n→∞ n→∞ 2n n→∞ 2n
ér

Sea r = lı́m an = lı́m bn . Luego tomando límite cuando n tiende a ∞ en la desigualdad

n→∞ n→∞
f (an ) f (bn ) < 0, se obtiene que ( f (r))2 ≤ 0. De allí que ( f (r))2 = 0, y en consecuencia
um

f (r) = 0, es decir, r es una raíz de f .

1
Por último sea cn = (an + bn ). Luego,
2
.N

1 1
|r − cn | ≤ |bn − an | ≤ n+1 (b0 − a0 ).
2 2
A

Ejemplo: Supongamos que se aplica el método de bisección para determinar una raíz r de
una función f en el intervalo [50, 63].

a) ¿Cuántos pasos se deberían realizar con el algoritmo de bisección si se desea obtener

una raíz en el intervalo [50, 63] con una precisión absoluta de 10−6 ?;

b) ¿Cuántos pasos se deberían realizar con el algoritmo de bisección si se desea obtener

una raíz en el intervalo [50, 63] con una precisión relativa de 10−6 ?

4
 Para responder a), y usando el Teorema anterior se tiene que
13
|r − cn | ≤ n+1
≤ 10−6 ,
2
y por lo tanto 13 106 ≤ 2n+1 . Luego tomando logaritmo natural a ambos lados, se tiene que

(n + 1) ln 2 ≥ ln(13) + 6 ln(10),

21
de donde se deduce que n ≥ 23.

20
Para el caso b), también se usa el teorema anterior y que 50 ≤ r ≤ 63, entonces

|r − cn | 1 13 13
≤ n+1 ≤ ≤ 10−6 .
|r| 2 |r| 50 2n+1

F
13 6
Luego 2n+1 ≥ 10 , y aplicando logaritmo natural a ambos miembros de la desigualdad

A
50
se deduce que n ≥ 17.

M
A
-F
o
ic
ér
um
.N
A

5
 Clase 5 - Solución de ecuaciones no lineales (2)

El problema

Recordemos el problema que comenzamos a estudiar en la clase anterior: dada f : R→R

no lineal, se desea hallar una solución r de la ecuación

21
f (x) = 0. (1)

20
Al igual que antes, estudiaremos otro método iterativo para resolver este problema, el cual
genere una sucesión de aproximaciones que se espera que converja a la solución buscada r.

F
A
Idea del método de Newton

M
Dado que en general la función f es no lineal, resolver la ecuación (1) es un problema
difícil. La idea del método de Newton consiste en reemplazar la resolución de este problema
A
difícil por la resolución de una sucesión de problemas fáciles, cuyas soluciones convergen
a la solución del problema difícil, bajo adecuadas hipótesis.
-F

Supongamos que r es una raíz de f y que x es una aproximación de r a una distancia h,

es decir, r = x + h o equivalentemente h = r − x. Supongamos también que f 00 existe y es
continua en un entorno de x que contiene a r, entonces haciendo el desarrollo de Taylor de f
centrado en x, tenemos que
o

0 = f (r) = f (x + h) = f (x) + f 0 (x)h + O(h2 ).

ic
ér

Cuando la aproximación x es próxima a r, el valor de h es pequeño y por lo tanto h2 es más

pequeño aún y así podríamos despreciar el término O(h2 )
um

0 = f (x) + h f 0 (x),

para despejar h
.N

f (x)
h=− .
f 0 (x)
A

Ahora bien, es claro que no conocemos exactamente h pues si tuviéramos x y h tendríamos la

solución r. Sin embargo, si la aproximación x está cerca de r, entonces la nueva aproximación
f (x)
(x + h) = x − 0 debería estar aún más cerca de r.
f (x)
Entonces, comenzando con una aproximación x0 de r, la iteración del método de Newton
consiste en calcular

f (xn )
xn+1 = xn − , n ≥ 0.
f 0 (xn )

1
 ê `Ě%႒ İ ʻǼεॷ Ž - \ é` Ě! ŠŽ

È Ü
ॷ ) ॷ͞ॷ ॷ 6ॷ ॷ::=B /ॷ ॷİॷ Ü/ॷ )
ॷ./ 'ॷ3ॷRŶ႒ 6, ॷÜ/ॷ3 ॷ
ॷ.
ॷ
͞ ӽǼKॷ
Ѭ İ JǼKॷ 6, ॷ͞ȸ JǼKॷ
Ѭ İȸ JǼK ¡ॷ €)—ûॷ )
ॷ Ü/
‫ॷڣ‬3 t /ॷॷ )
ॷ 
ॷ.—
ॷ 6, ॷ

ॷ 
ॷ:
ॷॷ İॷ  ॷ)
ॷÝ/ ॷ Rě႒ |ॷ/ॷÈ
 /ƻ  ॷ
) ॷ
ॷņ
ॷ/ॵµ Ü/Ôॷ 
ॷ
,
/ ॷ ॷ )
ॷ fŽՒµ =.
ॷ ॷ6ॷ:/ ॷ/
=ॷ <˅‫ ݩ‬6, ॷX/ Ü/ॷ)
=
ॷ 
ॷ@tÔ
t ॷ‫ك‬Üt
ॷ‫ݫ‬/
ņ
µॷ
Gráficamente, dado el)
ॷÏ punto (xn , f (xn ))' la
ŒB¡ॷ J|

ॷ—=
ॷ
idea consiste en aproximar el gráfico de la función
/7/ù˪
f por la recta tangenteБ

:tॷ
a f que/ॷpasa por (x n , f (xn/
/ ॷ )Üॷ =:)‫ॷھ‬ )).
=:=

Talrecta tangente
Ü/ûॷ 
ॷ está
½6,ॷ
Ü'ॷ dada‫ܫ‬/
ÜÔ/
ॷ ln (x) =
porÜ/ॷ
0
f (xn )(x − xn ) + f (xn6,). ॷVer
= Ü/ॷ:Ü/
Figura ॷ3=ॷ))ॷ
(1). )
ॷÈ
 /ॷÜ
= / İСɖިʩ¡ॷ | )ॷ6ॷƪ/ /ॷÜॷƜ/ॷ
/ॷÂÜ =
ॷ 'Ń*/˪ \/ॷ )ॷ
B:
ûॷ )
ॷ ):
ॷ3ॷ )
ॷ—=.
ॷॷ—)ॷ ) ॷ3=ॷ
ņ Ü/ॷ @ljú
Ü/ॷ.—
_ॷ )
ॷ
I
/
ॷ ˆ@ˠ Î˪ .
=
¡ॷ €) _ॷ/'ॷ3=ॷ  

/ ॷ& ॷÈ
/Ӱॷ

ωª“µӌ

21
!!.7!7
߂+úಏ By+ú© ಏ$࠘+ú©+úಏ‫ ڍ‬ఀ©ಏ

20
F
A
ɽ৆B O CMಏ <
Ě/
ƿ
/DD//
ƿ
Ŋ /Å  ƿ/ ƿ
M
A
Figura 1: Interpretación geométrica del método de Newton.
-F

Así, en lugar de buscar la raíz de f (problema difícil), se calcula la raíz de ln , es decir, se

busca xn+1 solución de ln (x) = 0. Es decir:
o

l(xn+1 ) = f 0 (xn )(xn+1 − xn ) + f (xn ) = 0,

ic

y por lo tanto
ér

f (xn )
xn+1 = xn − .
f 0 (xn )
um

Existen ejemplos donde esta idea podría fallar como se puede ver en la Figura (2). Mirando
el gráfico se puede deducir que la aproximación inicial x0 debe estar suficientemente cerca
de la raíz r para que el método converja.
.N

( Ɉ   Ȉ7 Ȉ
P
Ȉ
ȈF Ȉ 
Ȉ

Вêಏ
A

̓țӥ
 %"<   <
LƜ+  ƿ
ƿƿ
‰ƿ ƿ +êಏ

sĽ
Ǧ
k ॷ Ǧ  ॷLŪL
ॷ@ॷ@ Ǧ: ॷ )@ ॷ ॷ Ď )˪ 
˪ ๖ා႒ ྃ႒ Œ
˪=ॷ k@ ॷ k
ॷ
=@:kॷ3ॷİॷ)Ȣॷ
Figura ႒ :=
=&

2: Ejemplo donde elॷ k@:
ļॷ
método de Newton podría diverger.

¬¬ ¬àœ V‡Íµ{µà

2
Èॷ
ॷ )@$Ūॷ ȢРì'ƒ
ॷ k
ॷ
ē=ॷ ॷ È
 Ϊॷ Ǧ
k ţॷ Ăgॷ 
 3˪ 
ॷ Ǧ
@ॷ k
ॷ
I—@  
ॷ
 Algoritmo del método de Newton

Dados los siguientes datos de entrada y parámetros algorítmicos: la aproximación inicial x0 ,

el número máximo de iteraciones permitido M, δ la tolerancia para el error e (en la variable
x) y ε la tolerancia para los valores funcionales.
input x0 , M, δ , ε

21
v ← f (x0 )
output 0, x0 , v

20
if |v| < ε then STOP (1)
for k = 1, 2, . . . , M do

F
x1 ← x0 − v/ f 0 (x0 )
v ← f (x1 )

A
output k, x1 , v

M
if |x1 − x0 | < δ or |v| < ε then STOP (2)
x0 ← x1
A
end
-F

Observaciones:
1. En el algoritmo aparecen dos paradas (STOP). La primera detecta que el punto inicial
o

satisface la tolerancia del valor funcional. La segunda detecta que alguno de los criterios de
ic

parada considerados se cumple. Además el algoritmo podría parar por la cantidad máxima
de iteraciones permitida.
ér

2. Se podría usar otro criterio de parada basado en que el error relativo sea pequeño en lugar
de considerar el error absoluto.
um

3. El algoritmo requiere subprogramas o procedimientos externos que calculen f (x) y f 0 (x).

4. sería conveniente controlar en el programa que f 0 (x0 ) 6= 0 en cada iteración.
.N

Análisis de errores
A

El siguiente resultado muestra que bajo ciertas hipótesis la sucesión generada por el mé-
todo de Newton converge local y cuadráticamente, es decir que si se comienza con una
aproximación suficientemente próxima a la solución el método converge cuadráticamente.

Teorema 1. Si f 00 es continua en un entorno de una raíz r de f y si y si f 0 (r) 6= 0 entonces

existe δ > 0 tal que si el punto inicial x0 satisface |r − x0 | ≤ δ luego todos los puntos de la
sucesión {xn } generados por el algoritmo del método de Newton satisfacen que |r − xn | ≤ δ
para todo n, la sucesión {xn } converge a r y la convergencia es cuadrática, i.e., existen una

3
 constante c = c(δ ) y un natural N tal que

|r − xn+1 | ≤ c|r − xn |2 , para n ≥ N.

Demostración. Consideremos el error en la aproximación xn en la iteración n:

en = r − xn . (2)

21
El error en la iteración (n + 1) es:
 
f (xn ) f (xn ) f (xn )

20
en+1 = r − xn+1 = r − xn − 0 = r − xn + 0 = en + 0 ,
f (xn ) f (xn ) f (xn )
Luego
en f 0 (xn ) + f (xn )

F
en+1 = (3)
f 0 (xn )

A
Si escribimos el desarrollo de Taylor de f alrededor de xn tenemos que
1

para algún ξn entre xn y xn + h.

M
f (xn + h) = f (xn ) + f 0 (xn )h + f 00 (ξn )h2 ,
2
A
Tomando h = en = r − xn , se tiene que xn + h = r y luego
-F

1
0 = f (r) = f (xn ) + f 0 (xn )h + f 00 (ξn )h2 ,
2
para algún ξn entre xn y r. Luego
o

1
en f 0 (xn ) + f (xn ) = − f 00 (ξn )h2 . (4)
ic

2
ér

Usando (3) y (4) obtenemos que

1 f 00 (ξn ) 2
en+1 = − e . (5)
2 f 0 (xn ) n
um

Para acotar esta expresión definimos la siguiente función para δ > 0

máx | f 00 (x)|
.N

1 |x−r|≤δ
c(δ ) = .
2 mı́n | f 0 (x)|
|x−r|≤δ
A

Como f 0 y f 00 son funciones continuas entonces | f 0 | y | f 00 | alcanzan sus valores extremos en

un intervalo cerrado y acotado alrededor de r. Luego, dado δ > 0, para todo par x y ξ tal que
|ξ − r| ≤ δ y |x − r| ≤ δ existe una constante c = c(δ ) tal que
1 f 00 (ξ )
≤ c(δ ).
2 f 0 (x)

1 f 00 (r)
Ahora notemos que si δ → 0 entonces c(δ ) → , que es un número finito porque
2 f 0 (r)
f 0 (r) 6= 0 y por lo tanto δ c(δ ) → 0 cuando δ → 0. Entonces elegimos δ suficientemente
pequeño tal que ρ = δ c(δ ) < 1.

4
 Supongamos que la aproximación inicial x0 es tal que e0 = |x0 − r| ≤ δ , y como ξ0 está entre
x0 y r, entonces |ξ0 − r| ≤ δ . Luego por la definición de c(δ ) tenemos que
1 f 00 (ξ0 )
≤ c(δ ).
2 f 0 (x0 )

Finalmente por (5)

21
1 f 00 (ξ0 )
|x1 − r| = |e1 | = |e0 |2 ≤ c(δ )|e0 |2 = c(δ )|e0 ||e0 | ≤ c(δ )δ |e0 | = ρ|e0 | < |e0 | ≤ δ .
2 f 0 (x0 )

20
Esto dice que x1 está a una distancia de r menor que δ y además que si la sucesión converge,
lo hace cuadráticamente. Repitiendo este argumento obtenemos que
|e1 | ≤ ρ|e0 |

F
|e2 | ≤ ρ|e1 | ≤ ρ 2 |e0 |
|e3 | ≤ ρ|e2 | ≤ ρ 3 |e0 |

A
..
.
En general, tenemos que
M |en | ≤ ρ n |e0 |.
Como 0 < ρ < 1, entonces lı́m ρ n = 0 y por lo tanto lı́m en = 0 y así lı́m xn = r.
A
n→∞ n→∞ n→∞
-F

Observación: el método de Newton tiene convergencia cuadrática local para determinar

raíces simples, bajo adecuadas hipótesis.
o

El siguiente resultado que sólo enunciaremos aunque su demostración puede consultarse en

ic

la Bibliografía muestra que bajo hipótesis de convexidad en todo el dominio el método de

Newton converje independientemente del punto inicial.
ér

R R
Teorema 2. Si f 00 es continua en , f es creciente y convexa en y tiene una raíz, entonces
esa raíz es única y la iteración de Newton convergerá a esa raíz independientemente del
um

punto inicial x0 .

Ejercicio: hallar una función que cumpla las hipótesis del teorema. Analizar la utilidad prác-
tica del teorema.
.N

√
Aplicación del método de Newton: Sea R > 0 y x∗ = R, entonces x∗ es una raíz de x2 −R =
0. Si se aplica el método de Newton a f (x) = x2 − R, se obtiene la siguiente fórmula de
A

iteración:  
1 R
xn+1 = xn + .
2 xn
Por ejemplo si R = 17 y se comienza con x0 = 4, se obtienen
x1 = 4.12
x2 = 4.123106
x3 = 4.1231056256177
x4 = 4.123105625617660549821409856,
que tiene 28 dígitos correctos.

5
 Clase 6 - Solución de ecuaciones no lineales (3)

Método de la secante

Continuamos con el mismo problema de resolver una ecuación no lineal. Hasta ahora vimos
el método de bisección y el método de Newton. Recordemos que la iteración del método de

21
Newton está dada por:
f (xn )
xn+1 = xn − 0 para n ≥ 0. (1)
f (xn )

20
Si bien este método tiene convergencia cuadrática local, tiene como desventaja que requiere
la evaluación de la derivada de la función f en cada iteración. Uno de los métodos más
conocidos que evita esto es el método de la secante.

F
La idea del método de la secante consiste en reemplazar f 0 (xn ) en la iteración de Newton

A
(1) por una aproximación dada por el cociente incremental, dado por la pendiente de la recta
secante que pasa por los puntos (xn , f (xn )) y (xn + h, f (xn + h)):

M
f 0 (xn ) ≈ an =
f (xn + h) − f (xn )
h
,
A
para algún h suficientemente pequeño.
-F

Para evitar evaluar f en un punto adicional (xn + h) se elige h = xn−1 − xn , entonces:

f (xn + h) − f (xn ) f (xn−1 ) − f (xn ) f (xn ) − f (xn−1 )

an = = = .
h xn−1 − xn xn − xn−1
o
ic

Ver Figura 1.

-O   ȈT Ȉ

Ȉ
ȈFÿ Ȉ
Ȉ

ér
um
.N

 %"<   <
2M}tM„fJ¾fyM„€„MDf}y¾
}U¾ŒMJDy¾tMa}L¾
A

Figura 1: Interpretación gráfica del método secante.


ࢳ
ࢳ& ࢳࢳ ࢳ&J
%<ࢳ‚
&ࢳࢳࢳ&

ࢳ
 ¬ ˪‹
­ 3˪  ࢳࢳ
©ø
૵
 ࢄࢳf4ࢳ7 D i ࢳ6Mࢳ႒

Así, la iteración del método secante consiste en:


Ě - -ă,™T
Ě
4 \ 
Ě‫ݩݣ‬ 3C
Ě  @Ě
f (xn )
×
ࢳ&: xn+1
%ࢳ n − f (x )−
= x
ࢳࢳf4 ࢳf
:Ž
(xn−1
para n ≥ 1, ࢳVࢳX€ࢳ
%
ࢳ)&
ϲ_ࢳV&%ࢳ
ࢳ)
u
)
n
႒ %Ljࢳ%
_1ࢳ
x −x n n−1

©ø
૵ ©
૵ h†  Ž
ϣȧ—Ѭ
yå %ࢳ
v

es decir,  Ě ӌ
v

xn − xn−1
o
ࢳࢳ&̰ ̰%
ࢳࢳ%Š
1ࢳ
ࢳ&
ࢳ&%ࢳࢳ 
)ࢳА
ࢳ  ­˪
˪

©૵{
Ěʲ4Ě 
Ě

n+1 = xn − f (xn ) para n ≥ 1.
#ʲ
xЯࢳV%ࢳࢳ (2)
f (x ) − f (x )
n n−1
,  UĚ
ѬĚ
{Ě š @Ě
GJ ႁ႒ D i ࢳ .Ǎ
4Ě  b Ú
Ě
v
f֨

1
T 
ࢳ
ࢳ (  ࢳVࢳ, TĚ &
4&
ࢳĚ )ࢳ, ›  ĖĚࢳ  ࢳ ࢳ%ࢳ63ࢳ&
þ
 Algoritmo del método de la secante

Dados los siguientes datos de entrada y parámetros algorítmicos: a la penúltima aproxima-

ción de r, b la última aproximación de r, el número máximo de iteraciones permitido M, δ
la tolerancia para el error e (en la variable x) y ε la tolerancia para los valores funcionales.
input a, b, M, δ , ε

21
f a ← f (a)
f b ← f (b)

20
output 0, a, f a
output 1, b, f b

F
for k = 2, 3, . . . , M do
if | f a| < | f b| do

A
a ↔ b; f a ↔ f b
endif
s ← (b − a)/( f b − f a)
M
A
b←a
-F

fb ← fa
a ← a− fa∗s
f a ← f (a)
o

output k, a, f a
ic

if |b − a| < δ or | f a| < ε then STOP

ér

end do
um

Observaciones:
1. En el algoritmo los puntos a y b pueden intercambiarse para lograr que | f (b)| ≤ | f (a)|.
Así, para el par {xn−1 , xn } se satisface que | f (xn )| ≤ | f (xn−1 )|, y para el par siguiente
.N

{xn , xn+1 } se tiene que | f (xn+1 )| ≤ | f (xn )|. Esto garantiza que la sucesión {| f (xn )|} es no
creciente.
2. El algoritmo se detiene por el número máximo de iteraciones permitidas, por satisfacer la
A

tolarancia para los valores funcionales, o por la tolerancia en la diferencia de dos aproxima-
ciones sucesivas.
3. En cuanto al análisis de errores, es posible probar que:
en+1 ≈ ceαn = (c enα−1 ) e1n ,
√
donde α = (1 + 5)/2 = 1.618334... Como 1 < α < 2 se dice que el método de la secante
tiene convergencia superlineal. Además, por recurrencia
2
en+1 ≈ ceαn ≈ c1+α eαn−1
√ √
donde α 2 = (1 + 5)2 /4 = (3 + 5)/2 = 2.618334..., esto dice que dos iteraciones de
método de la secante es mejor que una iteración del método de Newton.

2
 Iteración de punto fijo

En esta sección veremos que es un punto fijo de una función dada, como encontrarlos numé-
ricamente y cuál es la conexión con el problema de determinar raíces de funciones.

21
Definición 1. un punto fijo de una función g es un número p, en el dominio de g, tal que
g(p) = p.

20
Por un lado, si p es una raíz de una función f , esto es, f (p) = 0, entonces es posible definir
diferentes funciones g con un punto fijo en p, por ejemplo: g(x) = x − f (x), o g(x) = x +
3 f (x).

F
Por otro lado, si g tiene un punto fijo en p, esto es, g(p) = p, entonces la función f (x) =

A
x − g(x) tiene una raíz en p.
Aunque estamos interesados en el problema de determinar soluciones de una ecuación no

M
lineal, o equivalentemente, encontrar raíces de funciones no lineales veremos que la forma de
la iteración de punto es muy fácil de estudiar y analizar. Además algunas opciones de punto
fijo dan origen a técnicas matemáticas y computacionales muy poderosas para determinar
A
raíces.
-F

Ejemplo 1: se busca determinar los posibles puntos fijos de la función g(x) = x2 − 2 en el

intervalo [−2, 3].
o

Para esto se plantea la ecuación g(x) = x, es decir, se debe resolver la ecuación cuadrática
x2 − 2 = x, o sea, x2 − x − 2 = 0. Luego es fácil verificar que x = −1 y x = 2 son puntos fijos
ic

(g(p) = p) pues
ér

g(−1) = (−1)2 − 2 = −1 y g(2) = (2)2 − 2 = 2.

Ver Figura 2.
um
.N
A

Figura 2: Puntos fijos de g.

A continuación veremos un resultado que da condiciones suficientes para la existencia y

unicidad del punto fijo.

3
 Teorema 1.

1. Si g ∈ C[a, b] (es decir, g es una función continua en el intervalo [a, b]) y g(x) ∈ [a, b]
para todo x ∈ [a, b] entonces existe p ∈ [a, b] tal que g(p) = p. (EXISTENCIA)

2. Si además existe g0 (x) para todo x ∈ (a, b) y existe una constante positiva k < 1 tal que
|g0 (x)| ≤ k para todo x ∈ (a, b), entonces el punto fijo en (a, b) es único. (Ver Figura 3).

21
(UNICIDAD)

20
F
A
M
Figura 3: Unicidad del punto fijo.
A
-F

Demostración.

1. Si g(a) = a o g(b) = b, el punto fijo está en uno de los extremos del intervalo y ya
o

estaría probado.
ic

Supongamos que esto no es cierto, entonces g(a) > a y g(b) < b. Sea h(x) = g(x) − x
una función definida en [a, b], que además es continua en [a, b] pues g y x lo son y resta
ér

de funciones continuas es una función continua. Además, por lo anterior, tenemos que

h(a) = g(a) − a > 0 y h(b) = g(b) − b < 0.

um

Entonces, por el Teorema del Valor Intermedio, existe p ∈ (a, b) tal que h(p) = 0, esto
es, g(p) = p y por lo tanto p es un punto fijo de g.
.N

2. Ahora supongamos que existen dos puntos fijos distintos p y q en [a, b], es decir, p, q ∈
[a, b], p 6= q tal que g(p) = p y g(q) = q. Por el Teorema del Valor Medio existe ξ
entre p y q, y por lo tanto en [a, b] tal que
A

g(p) − g(q) = g0 (ξ )(p − q).

Luego usando la hipótesis que |g0 (x)| ≤ k < 1 tenemos que

|p − q| = |g(p) − g(q)| = |g0 (ξ )||p − q| ≤ k|p − q| < |p − q|,

lo que es una contradicción que provino de suponer que habrían dos puntos fijos dis-
tintos en [a, b], y por lo tanto el punto fijo es único.

4
 Analicemos la existencia y unicidad de punto fijo en los siguientes ejemplos.
(x2 − 1) x2 1
Ejemplo 1: considerar g(x) = = − en el intervalo [−1, 1].
3 3 3
Es fácil ver g tiene un mínimo absoluto en x = 0 y g(0) = −1/3. Además tiene máximos
absolutos en x = −1 y x = 1 donde g(−1) = 0 y g(1) = 0. Además, claramente g es continua
2 2
en [−1, 1] y |g0 (x)| = x ≤ para todo x ∈ [−1, 1]. Por lo tanto existe un único punto fijo

21
3 3
p de g en el intervalo [−1, 1]. Para determinar p, planteamos

20
p2 − 1
g(p) = p, ⇒ = p, ⇒ p2 − 3p − 1 = 0,
3
√ √
(3 − 13) (3 13)

F
cuyas raíces son p1 = = −0,302776... y p2 = = 3,302776... El único
2 2
punto fijo es claramente p1 pues este punto está en [−1, 1] en cambio p2 no. Ver Figura 4.

A
M
A
-F
o
ic
ér
um

Figura 4: Puntos fijos del Ejemplo 1.

Notar en el gráfico que p2 también es un punto fijo de g en el intervalo [3, 4]. Sin embargo,
8
.N

g(4) = 5 y g0 (4) = > 1, por lo que no se satisfacen las hipótesis del teorema anterior en
3
el intervalo [3, 4]. Esto dice que tales hipótesis son condiciones suficientes para garantizar la
existencia y unicidad de un punto fijo, pero no son necesarias.
A

Ejemplo 2: considerar g(x) = 3−x = e−(ln 3)x en el intervalo [0, 1].

Su derivada es g0 (x) = −(ln 3)e−(ln 3)x = −(ln 3)3−x < 0, por lo tanto g es decreciente en el
intervalo [0, 1]. Entonces,
1
g(0) = 1 ≥ g(x) ≥ = g(1),
3
y por el teorema anterior, existe un punto fijo en [0, 1]. Por otro lado, g0 (0) = − ln 3 =
−1.0986... y por lo tanto no se puede usar el teorema anterior para garantizar unicidad pues
|g0 (x)| ≮ 1 en el intervalo (0, 1). Ver Figura 5.

5
 21
20
Figura 5: Puntos fijos del Ejemplo 1.

F
Idea del algoritmo de punto fijo

A
Para calcular aproximadamente el punto fijo de una función g primero se inicia con una
aproximación inicial p0 y calculando pn = g(pn−1 ) para n ≥ 1 se obtiene una sucesión de

M
aproximaciones {pn }. Si la función g es continua y la sucesión converge entonces lo hace a
un punto fijo p de g pues
A
p = lı́m pn = lı́m g(pn−1 ) = g( lı́m pn−1 ) = g(p).
n→∞ n→∞ n→∞
-F

Algoritmo del método de punto fijo

o

Dados los siguientes datos de entrada y parámetros algorítmicos: p0 una aproximación ini-
cial, el número máximo de iteraciones permitido M y δ la tolerancia para el error e (en la
ic

variable x)
input p0 , M, δ
ér

output 0, p0
um

i←1
while i ≤ M then do
p ← g(p0 )
.N

output i, p
if |p − p0 | < δ then STOP
A

i ← i+1
p0 ← p
end while

Observaciones:
1. El algoritmo es muy fácil de implementar.
2. los criterios de parada utilizados son la distancia entre dos iteraciones sucesivas y el nú-
mero máximo de iteraciones.

6
 Teorema 2.

Sea g ∈ C[a, b] tal que g(x) ∈ [a, b] para todo x ∈ [a, b]. Supongamos que existe g0 (x) para
todo x ∈ (a, b) y existe una constante positiva 0 < k < 1 tal que |g0 (x)| ≤ k para todo x ∈
(a, b), entonces para cualquier p0 ∈ [a, b] la sucesión definida por pn = g(pn−1 ), para n ≥ 1,

21
converge al único punto fijo p en (a, b).

Demostración. Por el teorema anterior, se sabe que bajo estas hipótesis existe un único punto

20
fijo p ∈ [a, b]. Como g(x) ∈ [a, b] para todo x ∈ [a, b], la sucesión de aproximaciones {pn }
está bien definida para todo n, es decir, pn ∈ [a, b] para todo n.
Para probar la convergencia se usará el Teorema del Valor Medio en lo siguiente

F
|pn − p| = |g(pn−1 ) − g(p)| = |g0 (ξn )||pn−1 − p| ≤ k|pn−1 − p|,

A
luego, por recurrencia, se tiene que

M
|pn − p| ≤ k|pn−1 − p| ≤ k2 |pn−2 − p| ≤ · · · ≤ kn |p0 − p|.
A
Como 0 < k < 1 entonces lı́m kn = 0, luego
n→∞
-F

lı́m |pn − p| ≤ lı́m kn |p0 − p| = |p0 − p| lı́m kn = 0,

n→∞ n→∞ n→∞

y por lo tanto, la sucesión {pn } converge al punto fijo p.

o

Corolario 1. Si g es una función que satisface las hipótesis del teorema anterior, se tienen
ic

las siguientes cotas de error

|pn − p| ≤ kn máx{p0 − a, b − p0 }
ér

kn
|pn − p| ≤ 1−k |p1 − p0 | para todo n≥1
um

Demostración. La demostración puede consultarse en el libro de Burden-Faires.

.N

Análisis de error en métodos de punto fijo

A

Supongamos que p es un punto fijo de una función g y que {pn } es la sucesión, que converge
a p, definida por pn+1 = g(pn ).
Sea pn una aproximación del punto fijo p, es decir pn = p + h, y consideremos el desarrollo
de Taylor de g centrado en p:
00 (p)
pn+1 = g(pn ) = g(p + h) = g(p) + g0 (p)(pn − p) + g 2 (pn − p)2 + . . .
(r−1) (r) (ξ )
(3)
(p)
+ g(r−1)! (pn − p)r−1 + g n
r! (pn − p)r ,

para algún ξn entre pn y p.

7
 Supongamos ahora que g0 (p) = g00 (p) = · · · = g(r−1) (p) = 0 pero g(r) (p) 6= 0, entonces

g(r) (ξn ) g(r) (ξn ) r

en+1 = pn+1 − p = g(pn ) − g(p) = (pn − p)r = en ,
r! r!
esto es,
g(r) (ξn ) r

21
en+1 = en ,
r!
y tomando límite se obtiene

20
|en+1 | |g(r) (p)|
lı́m = = C,
n→∞ |en |r r!

F
por lo tanto el método tiene orden de convergencia (al menos) r.
Conclusión: si las derivadas de la función de iteración de punto fijo se anulan en el punto

A
fijo p hasta el orden (r − 1) entonces el método tiene orden de convergencia (al menos) r.

M
Usando este resultado se obtienen tres corolarios interesantes que relacionan el método de
A
Newton con el método de punto fijo para una función f que satisface las hipótesis del teorema
de convergencia de Newton.
-F

Corolario 2. Si f es una función que tiene una raíz simple p, entonces el método de Newton
es un método de punto fijo y tiene orden de convergencia (al menos) 2.
o
ic

f (x)
Demostración. Sea g(x) = x − , la función de iteración del método de Newton. Es claro
f 0 (x)
ér

que si p es una solución de f (x) = 0, entonces p es un punto fijo de g pues

f (p)
um

g(p) = p − = p,
f 0 (p)

ya que f (p) = 0 y f 0 (p) 6= 0.

Ahora calculemos g0 (x) y luego evaluemos en p:
.N

( f 0 (x))2 − f 00 (x) f (x) f 00 (x) f (x) f 00 (x) f (x)

g0 (x) = 1 − = 1 − 1 + = ,
( f 0 (x))2 ( f 0 (x))2 ( f 0 (x))2
A

entonces
f 00 (p) f (p)
g0 (p) = = 0,
( f 0 (p))2
y por lo tanto el método tiene orden de convergencia (al menos) 2.

Corolario 3. Si p es una raíz de multiplicidad r ≥ 2 de f , entonces el método de Newton

tiene orden 1.

8
 f (x) 0 f (x) f 00 (x)
Demostración. Ya vimos que si g(x) = x − entonces g (x) = .
f 0 (x) ( f 0 (x))2
Ahora, supongamos que p es una raíz de multiplicidad r de f , esto es, f (x) = (x − p)r h(x),
con h una función tal que h(p) 6= 0 y r ≥ 2.
La derivada primera de f es

21
f 0 (x) = r(x − p)r−1 h(x) + (x − p)r h0 (x) = (x − p)r−1 rh(x) + (x − p)h0 (x) ,
 

y la derivada segunda de f es

20
f 00 (x) = r(r − 1)(x − p)r−2 h(x) + 2r(x − p)r−1 h0 (x) + (x − p)r h00 (x),
= (x − p)r−2 r(r − 1)h(x) + 2r(x − p)h0 (x) + (x − p)2 h00 (x) .
 

F
Luego

A
(x − p)r h(x)(x − p)r−2 r(r − 1)h(x) + 2r(x − p)h0 (x) + (x − p)2 h00 (x)
 
0
g (x) = ,

entonces M (x − p)2r−2 [rh(x) + (x − p)h0 (x)]2

h(p)r(r − 1)h(p) r − 1
A
g0 (p) = = 6= 0,
r2 (h(p))2 r
-F

pues r ≥ 2.

Por último, modificando levemente la función de iteración del método de Newton se puede
o

recuperar la convergencia cuadrática aún en casos de raíces múltiples.

ic

Corolario 4. Si p es una raíz de multiplicidad r ≥ 2 de f , entonces la siguiente modificación

ér

del método de Newton recupera la convergencia cuadrática:

f (xn ) f (x)
xn+1 = xn − r , esto es, g(x) = x − r .
um

f 0 (xn ) f 0 (x)

Demostración. Usando las expresiones de f , f 0 y f 00 obtenidas en el corolario anterior, se

tiene que
.N

( f 0 (x))2 − f 00 (x) f (x) f (x) f 00 (x)

g0 (x) = 1−r = 1−r+r 0
( f 0 (x))2 ( f (x))2
r r−2 r(r − 1)h(x) + 2r(x − p)h0 (x) + (x − p)2 h00 (x)
 
(x − p) h(x)(x − p)
A

= 1−r+r
(x − p)2r−2 [rh(x) + (x − p)h0 (x)]2
h(x) r(r − 1)h(x) + 2r(x − p)h0 (x) + (x − p)2 h00 (x)
 
= 1−r+r .
[rh(x) + (x − p)h0 (x)]2

Evaluando en el punto fijo p se obtiene

h(p)r(r − 1)h(p)
g0 (p) = 1 − r + r = 1 − r + (r − 1) = 0,
r2 (h(p))2

y por lo tanto el método de Newton modificado tiene convergencia (al menos) cuadrática.

9
 Clase 7 - Interpolación polinomial

El problema

Dada una tabla de (n + 1) puntos: (x0 , y0 ), (x1 , y1 ), . . . , (xn yn ) donde x0 , x1 , . . . , xn son distin-
tos, se desea determinar un polinomio p, con el menor grado posible, tal que

21
p(xi ) = yi para i = 0, . . . , n.
En este caso se dice que tal polinomio p interpola el conjunto de puntos dados. Ver Figura 1.

20
F
A
M
A
Figura 1: Interpolación polinomial.
-F

A continuación veremos un resultado de existencia y unicidad del polinomio interpolante.

Teorema 1. Dados x0 , x1 , . . . , xn números reales distintos con valores asociados y0 , y1 , . . . , yn
o

entonces existe un único polinomio pn de grado menor o igual a n tal que pn (xi ) = yi , para
i = 0, . . . , n.
ic

Demostración.
ér

Existencia: probaremos la existencia del polinomio interpolante p por inducción.

Para n = 0, es el caso obvio pues el polinomio constante p0 (x) = y0 satisface que tiene grado
um

menor o igual a 0 y que p0 (x0 ) = y0 .

Ahora supongamos, por hipótesis inductiva que se tiene un polinomio pk−1 de grado ≤ k − 1
con pk−1 (xi ) = yi para i = 0, . . . , k − 1. Vamos a construir el polinomio pk de grado ≤ k, tal
.N

que pk (xi ) = yi para i = 0, . . . , k, de la forma

pk (x) = pk−1 (x) + c(x − x0 )(x − x1 ) . . . (x − xk−1 ),
donde c es una constante a determinar. Notar que este polinomio tiene grado ≤ k. Además
A

pk interpola los primeros k puntos que interpola pk−1 pues

pk (xi ) = pk−1 (xi ) = yi , i = 0, . . . , k − 1.

El coeficiente c se determina usando la condición de interpolación pk (xk ) = yk , es decir:

pk (xk ) = pk−1 (xk ) + c(xk − x0 )(xk − x1 ) . . . (xk − xk−1 ) = yk ,
de donde se deduce que
yk − pk−1 (xk )
c= .
(xk − x0 )(xk − x1 ) . . . (xk − xk−1 )

1
 El coeficiente c está bien definido porque los números x0 , x1 , . . . xn son distintos y el denomi-
nador nunca se anula. Esto prueba la existencia del polinomio interpolante pk .
Unicidad: supongamos que existen dos polinomios interpolantes pn y qn de grado ≤ n, esto
es, pn (xi ) = yi y qn (xi ) = yi para i = 0, . . . , n.
Sea h = pn − qn . Claramente h es un polinomio de grado ≤ n. Además, h(xi ) = 0 para i =
0, . . . , n, por lo tanto h es un polinomio de grado ≤ n y tiene (n + 1) raíces reales. Luego, por

21
el teorema fundamental del Álgebra, h(x) = 0 para todo x y por lo tanto pn = qn .

20
Forma de Newton del polinomio interpolante

F
Si bien el polinomio interpolante es único, puede escribirse de diferentes formas. La forma
de Newton se obtiene inductivamente, como se hizo en la prueba de la existencia del teorema

A
anterior, agregando un nuevo término al polinomio interpolante de un grado menor.
Si n = 0: vimos que es suficiente definir el polinomio constante p0 (x) = c0 = y0 .
M
Si n = 1: dados los puntos (x0 , y0 ), (x1 , y1 ), se construye p1 tal que p1 (x) = c0 + c1 (x − x0 )
y p1 (x0 ) = c0 = y0 . Usando que p1 (x1 ) = y1 , entonces y1 = c0 + c1 (x1 − x0 ) y por lo tanto,
A
y1 − c0
c1 = .
x1 − x0
-F

Si n = k: tenemos que

pk (x) = pk−1 (x) + ck (x − x0 )(x − x1 ) . . . (x − xk−1 ),

o

y por recurrencia obtenemos que

ic

pk (x) = c0 + c1 (x − x0 ) + c2 (x − x0 )(x − x1 ) + · · · + ck (x − x0 ) . . . (x − xk−1 ).

ér

La forma de Newton compacta del polinomio interpolante resulta en

k i−1
um

pk (x) = ∑ ci ∏ (x − x j ).
i=0 j=0

m
Aquí se adopta la convención que ∏ (x − x j ) = 1 si m < 0.
.N

j=0

Para evaluar pk (x), una vez calculados los coeficientes ck , es conveniente usar el algoritmo
de Horner de multiplicación encajada.
A

Forma de Lagrange del polinomio interpolante

Veremos otra forma alternativa de expresar el polinomio interpolante pn , de grado ≤ n, aso-

ciado a los puntos (x0 , y0 ), (x1 , y1 ), . . . (xn , yn ) donde x0 , x1 , . . . , xn son distintos.
Primero se definen los polinomios básicos de Lagrange asociado a los puntos distintos
x0 , x1 , . . . , xn :
n (x − x )
j
li (x) = ∏ para i = 0, . . . , n.
(x
j=0 i − x j)
j6=i

2
 n (x − x j )
Así, por ejemplo, para i = 0, se tiene l0 (x) = ∏ .
j=1 (x0 − x j )

Es claro que el grado de li (x) es igual n para i = 0, . . . , n y que


1 si i = j
li (x j ) = δi j =
0 si i 6= j.

21
Ahora definimos la forma de Lagrange del polinomio interpolante por
n

20
pn (x) = ∑ yi li (x).
i=0

Es claro que pn es un polinomio de grado ≤ n y que pn (xi ) = yi , para i = 0, . . . , n.

F
n

A
Ejercicio: probar que ∑ li(x) = 1.
i=0

Sea p(x) =
n

i=0 M
∑ li(x). Es claro que p tiene grado ≤ n. Sea h(x) = p(x) − 1, un polinomio de
grado ≤ n. Además, h(xk ) = 0 para k = 0, . . . , n, es decir, h es polinomio de grado ≤ n y tiene
A
n + 1 raíces. Luego h(x) = 0 para todo x, y por lo tanto, p(x) = 1 para todo x.
-F

n
Esto también puede ser generalizado a ∑ ximli(x) = xm, si m ≤ n.
i=0
o

Error en el polinomio interpolante

ic

Ahora veremos un resultado sobre el error que se comete al reemplazar una función por un
polinomio que la interpola en algunos puntos dados. Primero veremos dos observaciones
ér

muy simples que serán útiles en el teorema siguiente.

Observación 1: si p es un polinomio de grado igual a 0 entonces p0 (x) ≡ 0;
um

si p es un polinomio de grado igual a 1 entonces p00 (x) ≡ 0;

si p es un polinomio de grado igual a 2 entonces p000 (x) ≡ 0; y en general,
si p es un polinomio de grado igual a n entonces p(n+1) (x) ≡ 0.
Observación 2: si f es una función continua en [a, b] y derivable en (a, b). Si además, f (a) =
.N

f (b) entonces existe α ∈ (a, b) tal que f 0 (α) = 0 (Teorema de Rolle). En particular, si f (a) =
f (b) = 0 entonces existe α ∈ (a, b) tal que f 0 (α) = 0. Más aún, si f (a) = f (b) = f (c) = 0
entonces existen α ∈ (a, b) y β ∈ (b, c) tal que f 0 (α) = f 0 (β ) = 0.
A

Teorema 2. Sea f una función en Cn+1 [a, b] y p un polinomio de grado ≤ n que interpola
a f en (n + 1) puntos distintos x0 , x1 , . . . , xn en [a, b]. Entonces para cada x ∈ [a, b] existe un
ξ = ξx ∈ (a, b) tal que
n
1 (n+1)
f (x) − p(x) = f (ξ ) ∏(x − xi ).
(n + 1)! i=0

Demostración. Si x = xi para i = 0, 1, . . . , n, la igualdad es trivialmente cierta, pues p inter-

pola a f en esos puntos, por hipótesis.

3
 Ahora supongamos que x 6= x0 , x1 , . . . , xn , y definimos
n
w(t) = ∏(t − xi) (polinomio en t)
i=0
f (x) − p(x)
c = (constante, pues x está fija)
w(x)
ϕ(t) = f (t) − p(t) − cw(t) (función en t),

21
la constante c está bien definida pues w(x) 6= 0 si x 6= xi para i = 0, . . . , n.
Notar que si t = x0 , x1 , . . . , xn entonces ϕ(t) = 0, pues p interpola a f en esos puntos y porque

20
w se anula allí.
Además, por las definiciones de w, c y ϕ es fácil ver que si t = x entonces ϕ(t) = 0.

F
Luego, ϕ(t) tiene (al menos) (n + 2) raíces: x0 , x1 , . . . , xn , x, y por el Teorema de Rolle (y la
Observación 2) tenemos que:

A
ϕ 0 (t) tiene (al menos) (n + 1) raíces en (a, b).

M
ϕ 00 (t) tiene (al menos) n raíces en (a, b).
ϕ 000 (t) tiene (al menos) (n − 1) raíces en (a, b).
A
..
.
-F

ϕ (n+1) (t) tiene (al menos) 1 raíz en (a, b). Sea ξ = ξx ∈ (a, b) tal raíz de ϕ (n+1) (t).
Luego,
0 = ϕ (n+1) (ξ ) = f (n+1) (ξ ) − p(n+1) (ξ ) − cw(n+1) (ξ ). (1)
o

n
Como ∏(t −xi ) = t n+1 + términos de orden menor, entonces w(n+1) (ξ ) = (n +1)!. Además,
ic

i=0
como p es un polinomio de grado menor o igual que n entonces, por la Observación 1,
p(n+1) (x) = 0 para todo x, en particular para x = ξ .
ér

Finalmente de (1) y de la definición de c,

um

f (x) − p(x)
0 = f (n+1) (ξ ) − c(n + 1)! = f (n+1) (ξ ) − (n + 1)!,
w(x)
y por lo tanto,
f (n+1) (ξ ) n
.N

f (x) − p(x) = (x − xi ).
(n + 1)! ∏ i=0
A

Ejemplo: dar una estimación del error que se comete al aproximar la función f (x) = sen x
por un polinomio interpolante de grado 9, que interpola a f en 10 puntos en el intervalo [0, 1].
Como el polinomio interpolante tiene grado n = 9 es claro que se requieren 10 puntos. Por el
teorema anterior, para poder acotar la expresión del teorema se necesita una cota de f (10) (ξ ).
Como f (x) = sen x, es fácil deducir que | f (10) (ξ )| ≤ 1.
Además, si bien no se conocen los puntos de interpolación, sabemos que todos ellos perte-
9
necen al intervalo [0, 1], por lo tanto ∏ |x − xi | ≤ 1.
i=0

4
 Entonces

f (n+1) (ξ ) n | f (n+1) (ξ )| n
1
| sen x − p(x)| = (x − x ) = ∏(x − xi) ≤ < 2.8 10−7 .
(n + 1)! ∏
i
i=0 (n + 1)! i=0 10!

21
Convergencia uniforme de los polinomios de interpolación

20
El teorema anterior da una expresión, para cada punto x, del error que se comete al interpolar
una función f por un polinomio interpolante p.
Ahora bien, si se construyen polinomios interpolantes pn utilizando cada vez más puntos,

F
o equivalentemente de grado cada vez más alto, sería natural esperar que estos polinomios
convergieran uniformemente a la función f en el intervalo [a, b]. Es decir, esperaríamos

A
k f − pn k∞ = máx | f (x) − pn (x)| → 0 si n → ∞.
a≤x≤b

M
Lamentablemente, para la mayoría de las funciones continuas no ocurre esto debido al com-
portamiento inestable de los polinomios de grado alto. Ver Figura 2.
A
-F
o
ic
ér
um
.N
A

Figura 2: Ejemplo de no convergencia uniforme de los polinomios interpolantes.

Como puede verse en la figura, a medida que se aumenta el grado del polinomio, y por lo
tanto la cantidad de puntos de interpolación, el gráfico del polinomio tiende a comportarse
muy diferente al gráfico de la función.

5
 Clase 8 - Interpolación polinomial (2)

Repaso:

- El problema: Dada una tabla de (n+1) puntos: (x0 , y0 ), (x1 , y1 ), . . . , (xn yn ), con x0 , x1 , . . . , xn
distintos, se desea determinar un polinomio p, con el menor grado posible, tal que

21
p(xi ) = yi para i = 0, . . . , n.

En este caso se dice que tal polinomio p interpola el conjunto de puntos dados.

20
- Existencia y unicidad del polinomio interpolante.
- Forma de Lagrange.
- Forma de Newton.

F
- Error en el polinomio interpolante.
- Convergencia de los polinomios de interpolación.

A
Diferencias divididas
M
Recordemos la forma de Newton del polinomio interpolante basado en los puntos distintos
A
x0 , x1 , . . . , xn :
-F

pn (x) = c0 + c1 (x − x0 ) + c2 (x − x0 )(x − x1 ) + · · · + cn (x − x0 ) . . . (x − xn−1 )

n i−1
= ∑ ci ∏ (x − x j ).
i=0 j=0
o

Notemos que el coeficiente c0 se obtiene simplemente de c0 = y0 (o de f (x0 )). El coeficiente

ic

c1 se calcula con y1 (o de f (x1 )), despejando en

y1 − y0
ér

y1 = p1 (x1 ) = c0 + c1 (x1 − x0 ), esto es c1 = ,

x1 − x0
este coeficiente depende de x0 , x1 , y0 , y1 . En general, para calcular el coeficiente ck se re-
um

quieren conocer x0 , . . . , xk , y0 , . . . , yk , o si estamos interpolando a una función f se requieren:

x0 , . . . , xk , f (x0 ), . . . , f (xk ). Este coeficiente se denota por

ck = f [x0 , x1 , . . . , xk ],
.N

para k = 0, . . . , n y se denomina diferencias divididas.

Ahora la forma de Newton compacta del polinomio interpolante resulta
A

k i−1
pk (x) = ∑ f [x0 , x1 , . . . , xi ] ∏ (x − x j ),
i=0 j=0

donde si k = 0:
f [x0 ] = f (x0 ),
si k = 1:
f (x1 ) − f (x0 )
f [x0 , x1 ] = .
x1 − x0
Veremos a continuación un resultado general para determinar las diferencias divididas aso-
ciadas a un polinomio que interpola una función f .

1
 Teorema 1. Dados x0 , x1 , . . . , xn números reales distintos, las diferencias divididas satisfacen
la siguiente ecuación

f [x1 , x2 , . . . , xn ] − f [x0 , x1 , . . . , xn−1 ]

f [x0 , x1 , . . . , xn ] = .
xn − x0

21
Demostración.
Sea pn−1 el polinomio de grado ≤ n − 1 que interpola a f en los n puntos x0 , x1 , . . . , xn−1 .

20
Sea q el polinomio de grado ≤ n − 1 que interpola a f en los n puntos x1 , x2 , . . . , xn .
Sea pn el polinomio de grado ≤ n que interpola a f en los n + 1 puntos x0 , x1 , . . . , xn .

F
Se afirma que
(x − xn )

A
pn (x) = q(x) + [q(x) − pn−1 (x)]. (1)
(xn − x0 )

M
Es claro que a ambos lados de la igualdad se tienen polinomios de grado ≤ n. Además,
para i = 0,
A
(x0 − xn )
pn (x0 ) = f (x0 ) y q(x0 ) + [q(x0 ) − pn−1 (x0 )] = pn−1 (x0 ) = f (x0 ).
-F

(xn − x0 )

Para i = 1, . . . , n − 1,
o

(xi − xn )
pn (xi ) = f (xi ) y q(xi ) + [q(xi ) − pn−1 (xi )] = f (xi ),
ic

(xn − x0 )
ér

pues q(xi ) = pn−1 (xi ) = f (xi ) para i = 1, . . . , n − 1.

Para i = n,
um

(xn − xn )
pn (xn ) = f (xn ) y q(xn ) + [q(xn ) − pn−1 (xn )] = q(xn ) = f (xn ).
(xn − x0 )
.N

Por lo tanto, a ambos lados lados de (1) se tiene un polinomio de grado ≤ n y ambos interpo-
lan a f en los mismos n + 1 puntos distintos. Luego por unicidad del polinomio interpolante,
la igualdad (1) es cierta. Como ambos polinomios son iguales, los coeficientes de cada po-
tencia de x deben coincidir. En particular considerando el coeficiente de xn a ambos lados de
A

(1) obtenemos

1
f [x0 , x1 , . . . , xn ] = ( f [x1 , x2 , . . . , xn ] − f [x0 , x1 , . . . , xn−1 ])
xn − x0
f [x1 , x2 , . . . , xn ] − f [x0 , x1 , . . . , xn−1 ]
= .
xn − x0

2
 Luego

f [x0 ] = f (x0 )

f [x1 ] − f [x0 ]
f [x0 , x1 ] =

21
x1 − x0
f [x1 , x2 ] − f [x0 , x1 ]
f [x0 , x1 , x2 ] =

20
x2 − x0
.. ..
. = .

F
f [xi+1 , xi+2 , . . . , xi+ j ] − f [xi , xi+1 , . . . , xi+ j−1 ]
f [xi , xi+1 , . . . , xi+ j ] = (2)
xi+ j − xi

A
Tabla de diferencias divididas
M
A
Dados 4 puntos distintos (no necesariamente ordenados) se puede construir la tabla de dife-
rencias divididas de la siguiente manera:
-F

x0 f [x0 ] f [x0 , x1 ] f [x0 , x1 , x2 ] f [x0 , x1 , x2 , x3 ]

x1 f [x1 ] f [x1 , x2 ] f [x1 , x2 , x3 ]
o

x2 f [x2 ] f [x2 , x3 ]
ic

x3 f [x3 ]
ér

Ejemplo: Dados los siguientes valores:

um

x 3 1 5 6
f (x) 1 −3 2 4
.N

a) obtener la tabla de diferencias divididas

3 1 2 −3/8 7/40
A

1 −3 5/4 3/20
5 2 2
6 4

b) obtener el polinomio interpolante en la forma de Newton:

3 7
p(x) = 1 + 2(x − 3) − (x − 3)(x − 1) + (x − 3)(x − 1)(x − 5).
8 40

3
 Algoritmo para calcular los coeficientes de la tabla de diferencias divididas

Para obtener algorítmicamente los coeficientes de la tabla de diferencias divididas se puede

pensar la misma como en un arreglo matricial de la siguiente forma:
x0 c00 c01 c02 c03 . . . c0,n−1 c0,n
x1 c10 c11 c12 c13 . . . c1,n−1

21
.
x2 c20 c21 c22 c23 . .
.. .. .. ..

20
.
. . . . ..
.. .. .. .
. . . ..

F
xn−1 cn−1,0 cn−1,1
xn cn,0

A
donde ci j = f [xi , xi+1 , . . . , xi+ j ].

M
Dados los datos de entrada x0 , x1 , . . . , xn y f (x0 ), f (x1 ), . . . , f (xn ) se pueden calcular estos
coeficientes a partir de la fórmula (2) de la siguiente manera:
A
for j = 1 to n do
-F

for i = 0 to n − j do
ci, j ← (ci+1, j−1 − ci, j−1 )/(xi+ j − xi )
end do
o

end do
ic

Este algoritmo puede ser optimizado de modo de almacenar estos coeficientes en un vector
en vez de una matriz para ahorrar espacio de almacenamiento.
ér
um

Propiedades de las diferencias divididas

Veremos a continuación algunos resultados interesantes acerca de las diferencias divididas.

.N

El primero es un resultado sobre la permutación de los puntos de interpolación, el segundo

sobre el error de interpolación y el tercero y último relaciona las diferencias divididas con
las derivadas de orden superior.
A

Teorema 2. Sean x0 , x1 , . . . , xn números reales distintos y z0 , z1 , . . . , zn un reordenamiento de

x0 , x1 , . . . , xn . Entonces f [z0 , z1 , . . . , zn ] = f [x0 , x1 , . . . , xn ].

Demostración. Primero es importante notar que el polinomio interpolante de grado ≤ n ba-

sado en los puntos x0 , x1 , . . . , xn coincide con el polinomio interpolante de grado ≤ n basado
en los puntos z0 , z1 , . . . , zn , por unicidad del polinomio interpolante.
Luego, el coeficiente de xn en el polinomio de grado ≤ n que interpola a f en z0 , z1 , . . . , zn
es f [z0 , z1 , . . . , zn ] y el coeficiente de xn en el polinomio de grado ≤ n que interpola a f en
x0 , x1 , . . . , xn es f [x0 , x1 , . . . , xn ], y deben ser iguales pues estos dos polinomios son iguales.

4
 Teorema 3. Sea p el polinomio de grado ≤ n que interpola a f en los n + 1 nodos distintos
x0 , x1 , . . . , xn . Si t es un número real distinto de los nodos, entonces
n
f (t) − p(t) = f [x0 , x1 , . . . , xn ,t] ∏ (t − x j ).

21
j=0

Demostración. Sea q el polinomio de grado ≤ n + 1 que interpola a f en x0 , x1 , . . . , xn ,t. Por

20
la manera en que se genera el polinomio interpolante en la forma de Newton se sabe que
n
q(x) = p(x) + f [x0 , x1 , . . . , xn ,t] ∏ (x − x j ).

F
j=0

A
Como q interpola a f en el punto t, se tiene que q(t) = f (t), y entonces:

M n
f (t) = p(t) + f [x0 , x1 , . . . , xn ,t] ∏ (t − x j ),
j=0
A
por lo tanto,
n
-F

f (t) − p(t) = f [x0 , x1 , . . . , xn ,t] ∏ (t − x j ).

j=0
o

Teorema 4. Si f es una función n veces continuamente diferenciable en [a, b] y x0 , x1 , . . . , xn

ic

son n + 1 nodos distintos en [a, b], entonces existe un punto ξ ∈ (a, b) tal que

1 (n)
ér

f [x0 , x1 , . . . , xn ] = f (ξ ).
n!
um

Demostración. Sea p el polinomio de grado ≤ n − 1 que interpola a f en x0 , x1 , . . . , xn−1 . Por

el teorema del error en el polinomio interpolante de la clase anterior, aplicado a x = xn , se
sabe que
n−1
1
.N

f (xn ) − p(xn ) = f (n) (ξ ) ∏ (xn − x j ).

n! j=0

Ahora, por el Teorema 3 se obtiene

A

n−1
f (xn ) − p(xn ) = f [x0 , x1 , . . . , xn ] ∏ (xn − x j ),
j=0

y por lo tanto
1 (n)
f [x0 , x1 , . . . , xn ] = f (ξ ).
n!

5
 Interpolación de Hermite

Comenzaremos analizando a las diferencias divididas como función de sus argumentos. Has-
ta ahora hemos definido a las diferencias divididas para puntos distintos x0 , x1 , . . . , xn . Con-
sideremos ahora el caso simple de 2 puntos x0 y x:

f [x] − f [x0 ] f (x) − f (x0 )

21
f [x0 , x] = = .
x − x0 x − x0
Ahora tomando límite para x que tiende a x0 se tiene que

20
f (x) − f (x0 )
lı́m f [x0 , x] = lı́m = f 0 (x0 ).
x→x0 x→x0 x − x0

F
Luego es posible extender la definición de diferencias divididas para números repetidos de
la siguiente manera

A
f (x) − f (x0 )
= f 0 (x0 ).
x→x0
M
f [x0 , x0 ] = lı́m f [x0 , x] = lı́m
x→x0 x − x0

Con esta generalización se pueden construir un polinomio interpolante de grado 3, con sólo
A
2 puntos de interpolación y agregando 2 condiciones de interpolación de la derivada en esos
-F

mismos puntos, esto es:

p(x0 ) = f (x0 ), p(x1 ) = f (x1 ),

p0 (x0 ) = f 0 (x0 ), p0 (x1 ) = f 0 (x1 ).
o
ic

Así se obtiene la siguiente tabla de diferencia dividida

x0 f [x0 ] f 0 (x0 ) f [x0 , x0 , x1 ] f [x0 , x0 , x1 , x1 ]

ér

x0 f [x0 ] f [x0 , x1 ] f [x0 , x1 , x1 ]

um

x1 f [x1 ] f 0 (x1 )
x1 f [x1 ]
.N

Ahora, el polinomio interpolante basado en esta tabla está dado por

p(x) = f [x0 ] + f 0 (x0 )(x − x0 ) + f [x0 , x0 , x1 ](x − x0 )2 + f [x0 , x0 , x1 , x1 ](x − x0 )2 (x − x1 ).

A

En general, el polinomio interpolante que usa las derivadas en un punto se llama forma de
Hermite.
El siguiente resultado es una generalización del Teorema 4.

Teorema 5. Se f una función definida en [a, b], n veces continuamente diferenciable en [a, b].
Sean x0 , x1 , . . . , xn ∈ [a, b] puntos distintos y z ∈ (a, b). Entonces

f (n) (z)
lı́m f [x0 , x1 , . . . , xn ] = .
(x0 ,x1 ,...,xn )→(z,z,...,z) n!

6
 Demostración. Por el Teorema 4 se sabe que existe ξ ∈ (a, b) tal que

1 (n)
f [x0 , x1 , . . . , xn ] = f (ξ ).
n!

Como x0 → z, x1 → z, . . . , xn → z entonces ξ → z.
Como f (n) es continua entonces

21
1 (n) f (n) (z)
lı́m f [x0 , x1 , . . . , xn ] = lı́m f (ξ ) = .

20
(x0 ,x1 ,...,xn )→(z,z,...,z) ξ →z n! n!

F
El siguiente corolario es una consecuencia directa del teorema anterior.

A
Corolario 1. Si f es n veces continuamente diferenciable en un entorno del punto x0 , enton-
ces
f (n) (x0 )

M
f [x0 , x0 , . . . , x0 ] =
n!
.
A
Demostración. Trivial, basta tomar z = x0 .
-F

Ejemplo: determinar el polinomio de grado 4 que interpola los siguientes datos:

p(1) = 2, p0 (1) = 3,
o

p(2) = 6, p0 (2) = 7, p00 (2) = 8.

ic

La tabla de diferencias divididas resulta en

ér

1 2 3 1 2 −1
1 2 4 3 1
um

2 6 7 4
2 6 7
2 6
.N

y el polinomio interpolante de Hermite es

p(x) = 2 + 3(x − 1) + (x − 1)2 + 2(x − 1)2 (x − 2) − (x − 1)2 (x − 2)2 .

A

7
 Clase 9 - Interpolación polinomial (3)

Repaso

El problema: Dada una tabla de (n+1) puntos: (x0 , y0 ), (x1 , y1 ), . . . , (xn yn ) donde x0 , x1 , . . . , xn
son distintos, se desea determinar un polinomio p, con el menor grado posible, tal que:

21
p(xi ) = yi para i = 0, . . . , n.

20
En este caso se dice que tal polinomio p interpola el conjunto de puntos dados.
- Existencia y unicidad del polinomio interpolante.
- Forma de Lagrange.
- Forma de Newton.

F
- Error en el polinomio interpolante.
- Convergencia de los polinomios de interpolación.

A
- Diferencias divididas.
- Polinomios de Hermite.

M
A
Splines
-F

Antes de introducir el concepto de splines vamos a considerar el caso simple de interpolación

lineal que será muy útil en lo que sigue.
Sea f una función definida en el intervalo [x0 , x1 ] 2 veces continuamente diferenciable. El
o

polinomio de grado ≤ 1 que interpola a f en los puntos x0 , x1 es:

ic

p(x) = f [x0 ] + f [x0 , x1 ](x − x0 ),

ér

y el error está dado por

f 00 (ξ )
um

e(x) = f [x0 , x1 , x](x − x0 )(x − x1 ) = (x − x0 )(x − x1 ),

2!
para x, ξ ∈ (x0 , x1 ).
Sea M > 0 una constante tal que | f 00 (x)| ≤ M para todo x ∈ [x0 , x1 ].
.N

Sea ϕ(x) = (x − x0 )(x − x1 ), una función cuadrática, cuyo gráfico es una parábola con las
ramas hacia arriba, sus raíces son x0 y x1 y su mínimo se alcanza en xm = (x0 + x1 )/2.
A

Por lo tanto
  
x0 + x1 x0 + x1
|ϕ(x)| ≤ |(xm − x0 )(xm − x1 )| = − x0 − x1
2 2
(x1 − x0 ) (x0 − x1 ) |x1 − x0 |2
= = .
2 2 4

Por lo tanto,
M
|e(x)| ≤ |x1 − x0 |2 . (1)
8

1
 Supongamos que se desea interpolar una función f por un polinomio interpolante pn . Usar
pocos puntos de interpolación podría generar un polinomio que no aproxime bien a la fun-
ción. Por otro lado, como se comentó anteriormente, y contrariamente a lo que podría espe-
rarse, aumentar la cantidad de puntos de interpolación no mejora la convergencia uniforme
del polinomio interpolante pn a la función f . Esto es conocido como fenómeno de Runge.
Una idea que trata de conciliar estos conceptos opuestos es aplicar interpolación con polino-

21
mios de grado bajo por subintervalos. Esto es conocido como interpolación polinomial por
partes o interpolación segmentaria o simplemente splines.

20
Definición 1. Una función spline está formada por polinomios definidos en subintervalos,
los cuales se unen entre sí obedeciendo ciertas condiciones de continuidad.
Más formalmente, dados n + 1 puntos tales que x0 < x1 < . . . < xn , que denominaremos
nodos, y un entero k ≥ 0, un spline de grado k es una función S definida en [x0 , xn ] que

F
satisface:

A
- S es un polinomio de grado ≤ k en cada subintervalo [xi , xi+1 ), para i = 0, . . . , n − 1;
- las derivadas S(i) son continuas en [x0 , xn ], para i = 0, . . . , k − 1.

M
Veremos con un poco más de detalles los splines lineales y cúbicos, esto es, de grado 1 y 3.
A
-F

Splines lineales

Dados los n + 1 nodos tales que x0 < x1 < . . . < xn , un spline lineal (k = 1) es una función S
definida en [x0 , xn ] que satisface:
o

- S es un polinomio de grado ≤ 1 (recta) en cada subintervalo [xi , xi+1 ), para i = 0, . . . , n − 1;

ic

- la función S es continua en [x0 , xn ].

Es decir,
ér



 S0 (x) = a0 x + b0 , x ∈ [x0 , x1 )
 S1 (x) = a1 x + b1 , x ∈ [x1 , x2 )

um

S(x) = .. .. ,


 . .
Sn−1 (x) = an−1 x + bn−1 , x ∈ [xn−1 , xn ]


donde los 2n coeficientes ai , bi , para i = 0, . . . , n − 1 son las incógnitas a ser determinadas.

.N

Para eso, se deben tener 2n condiciones.

@8Ɉ 8÷¾# Ȉ ‘$Ȉ ¾¾B$ƟO Ȉm O¹'O=Ȉ
Notar que la segunda condición significa que los(J polinomios de grado ≤ 1 se pegan bien en
. .( J, x‫ݥ‬႒ n−1 . Las (n؊؋
los (n − 1) nodos interiores x"5D1 , . ‫ݤݣ‬ ،ಏ
ʖ႒ + 1) condiciones faltantes corresponden a las
H႒
A

(n + 1) condiciones de interpolación ' S(x $' i ) = f H႒(xi )

ј႒ para i = 0, . . . , n. Ver Figura 1.
΁ ú ‫ݩ‬
' ; 3J H႒ ʖ႒
;J ;J
H႒ H႒
TcX zrOÛ7!.Û "D
ÄâXâñĕ"
௜  ͫ


௜@ͫ (DJ

S]X zrOÛ7"1Û
Äââĕ
௜Ƣ


௜ #Ƣ

' B Bߠ ×řϼߠ F.ߠ Ǥ
q
ߠ Bߠ 
ߠ
B
ߠ
ߠ Ʌ B
Gߠ Ž ߠ ׏B.B .B޿.&ߠ 
ߠ ×řBߠ 
ߠ 
ߠ

;
 Ͻߠ 4ΝăяBߠ
ߠFigura
™ğĨ B# ߠ ώ -1:m)yGráfico de spline
†i $Ń0 ߠ
ߠ
; ™Ġ€  Ĩ

 Bߠ lineal-(k
6 )E= 1)
Bߠ#ߠ
B
ߠ 1
„UĜ 7Ńm7y :
ߠ' B Bߠ Ń ߠB B  &ߠřB ߠߠ
ߠ


ߠBDZ ߠ
ߠ Eߠߠ
ߠ B
Ƒߠߠ &ߠ  
  Ć $ Ń%  M 
  i $  Ń ¨ॷ 
ζ
ߠߠ
0 
ߠ
ߠ
q Ð

#ߠ B
ߠ *$Ń' ͈ ߠ
2
{„ƨ|`Ǎ D;  QϏ౱ Ӑ ֶ ಏ3;  Q¶ ʳ ‫ ݩ‬႒
+૵

˪ % Êॷ 
˪19 ْ Êॷ 
˪
 Dado un i fijo, se pueden determinar los coeficientes ai , bi , para i = 0, . . . , n−1 de la siguiente
manera:
ai xi + bi = Si (xi ) = f (xi )
ai xi+1 + bi = lı́m Si (x) = Si+1 (xi+1 ) = f (xi+1 )
x→xi+1

Restando la segunda ecuación menos la primera obtenemos ai (xi+1 − xi ) = f (xi+1 ) − f (xi ),

21
y por lo tanto

f (xi+1 ) − f (xi )

20
ai = , bi = f (xi ) − ai xi .
(xi+1 − xi )

Observación: supongamos que f es 2 veces continuamente diferenciable en [a, b] y xk =

a + kh, k = 0, . . . , n, con h = (b − a)/n.

F
Si S es un spline lineal, en cada intervalo [xk , xk+1 ] se tiene un polinomio de grado ≤ 1.

A
Entonces el error de interpolación para cada x ∈ [a, b] está dado por:
M 2

M |e(x)| ≤

donde | f 00 (x)| ≤ M para todo x ∈ [a, b] = [x0 , xn ].

8
h ,
A
-F

Splines cúbicos

Dados los n + 1 nodos tales que x0 < x1 < . . . < xn , un spline cúbico (k = 3) es una función
S definida en [x0 , xn ] que satisface:
o

S es un polinomio de grado ≤ 3 en cada subintervalo [xi , xi+1 ), para i = 0, . . . , n − 1;

ic

las funciones S, S0 y S00 son continuas en [x0 , xn ].

ér

Es decir,
S0 (x) = a0 x3 + b0 x2 + c0 x + d0 ,

x ∈ [x0 , x1 )
um


 S1 (x) = a1 x3 + b1 x2 + c1 x + d1 ,

x ∈ [x1 , x2 )

S(x) = .. .. ,


 . .
Sn−1 (x) = an−1 x3 + bn−1 x2 + cn−1 x + dn−1 , x ∈ [xn−1 , xn ]

.N

donde los 4n coeficientes ai , bi , ci , di , para i = 0, . . . , n − 1 son las incógnitas a ser determina-

das. Para eso, se deben tener 4n condiciones.
A

S(xi ) = f (xi ), i = 0, . . . , n ((n+1) condiciones de interpolación)

Si (xi+1 ) = lı́m Si (x) = Si+1 (xi+1 ), i = 0, . . . , n − 2 ((n-1) condiciones de continuidad de S)
x→xi+1
Si0 (xi+1 ) = lı́m Si0 (x) = Si+1
0
(xi+1 ), i = 0, . . . , n − 2 ((n-1) condiciones de continuidad de S0 )
x→xi+1
Si00 (xi+1 ) = lı́m Si00 (x) = Si+1
00
(xi+1 ), i = 0, . . . , n − 2 ((n-1) condiciones de continuidad de S00 )
x→xi+1

Esto da un total de (n + 1) + 3(n − 1) = 4n − 2 condiciones. Para poder determinar una única

solución se deben imponer dos condiciones adicionales:
S00 (x0 ) = S000 (x0 ) = 0 y S00 (xn ) = Sn−1
00
(xn ) = 0.

3
 En este caso, se denomina spline con condiciones naturales o simplemente spline natural.
Otras veces se suele utilizar

S0 (x0 ) = S00 (x0 ) = f 0 (x0 ) y S0 (xn ) = Sn−1

0
(xn ) = f 0 (xn ).

En este caso, se denomina spline con condiciones correctas.

21
Estas condiciones suelen estar asociadas a características del problema y son indicadas en el
problema o proporcionadas por quien presenta el problema.

20
F
A
M
A
-F
o
ic
ér
um
.N
A

4
 Clase 10 - Aproximación de funciones

El estudio de aproximación de funciones o, mejor dicho, de la teoría de aproximación com-

prende dos tipos de problemas: i) la búsqueda de parámetros óptimos de un modelo funcional
propuesto que represente un conjunto de datos dados; y ii) dada una función de manera ex-
plícita, se desea encontrar un tipo o modelo más sencillo para representarla, por ejemplo un

21
polinomio, que permita estimar valores funcionales de una manera más simple.

20
Aproximación discreta por cuadrados mínimos

Consideremos el problema de estimar una función desconocida a través de un conjunto de

F
datos experimentales: (xi , yi ), i = 1, . . . , m.
Por ejemplo, consideremos los siguientes datos

A
xi
yi
1
1.3
2
3.5
3
4.2M 4
5.0
5
7.0
6
8.8 10.1
7 8
12.5
9
13.0
10
15.6
A
Podemos representarlos gráficamente
-F
o
ic
ér
um
.N

Figura 1: Datos medidos

A

Si se construyera un polinomio interpolante usando los m puntos se obtendría un polinomio

de grado alto y ya se mencionó que esto no es conveniente. Para el caso particular de los
datos del ejemplo dado donde se tienen 10 puntos se obtendría un polinomio interpolante de
grado 9. Ver Figura (2).
Sin embargo, no es razonable exigir que la función pase exactamente por todos esos puntos.
La Figura (1) sugiere que la relación entre x e y es lineal (polinomio de grado 1). El motivo
que no exista una recta que se ajuste a estos datos, es decir que pase por todos los puntos, es
que estos datos tienen errores.

1
 21
20
F
A
Figura 2: Polinomio de grado 9

M
Sería mejor pensar en un modelo lineal y encontrar la mejor recta (en cierto sentido), aún
cuando no coincida con los datos en ningún punto. El problema ahora será determinar cuál
A
es la mejor recta
-F

y(x) = a1 x + a0
que mejor ajusta a esos datos, esto es, determinar los coeficientes a1 y a0 óptimos.
Una idea posible es determinar a0 y a1 tales que minimicen la función
o

E∞ (a0 , a1 ) = máx |yi − (a1 xi + a0 )|.

1≤i≤10
ic

Esto es conocido como el problema minimax.

ér

Otra alternativa es determinar a0 y a1 tales que minimicen la función

10
um

E1 (a0 , a1 ) = ∑ |yi − (a1 xi + a0 )|.

i=1

Esa función es conocida como desviación absoluta.

.N

Por último, el método de cuadrados mínimos para ajustar a una recta con m datos consiste
en determinar a0 y a1 tales que minimicen la función
m
A

E(a0 , a1 ) = E2 (a0 , a1 ) = ∑ [yi − (a1 xi + a0 )]2

i=1

con respecto a las variables a0 y a1 . En el ejemplo anterior m = 10.

Una condición necesaria para tener un mínimo es que las derivadas parciales de E con res-
pecto a a0 y a1 deben ser cero, esto es,

∂ ∂ m 2
m
∂ a0 E(a0 , a1 ) = ∑ [yi − (a 1 xi + a0 )] = 2 ∑ (yi − a1xi − a0)(−1) = 0,
∂ a0 i=1 i=1
∂ ∂ m 2
m
∂ a1 E(a 0 , a 1 ) = ∑ [yi − (a 1 xi + a0 )] = 2 ∑ (yi − a1 xi − a0 )(−xi ) = 0.
∂ a1 i=1 i=1

2
 Luego
m m m
(−2) ∑ (yi − a1 xi − a0 ) = 0 ⇒ ∑ yi − a1 ∑ xi − m a0 = 0,
i=1 i=1 i=1
m m m m
(−2) ∑ (yi − a1 xi − a0 )xi = 0 ⇒ ∑ xiyi − a1 ∑ xi2 − a0 ∑ xi = 0.
i=1 i=1 i=1 i=1

21
Reordenando estas ecuaciones se puede obtener el siguiente sistema lineal de dos ecuaciones
con las 2 incógnitas a0 y a1 :

20
 m m



 a 0 m + a1∑ i x = ∑ yi
i=1 i=1
m m m
2
a x + a1 ∑ xi = ∑ xiyi.

 0∑ i



F
i=1 i=1 i=1

A
Estas ecuaciones son conocidas como ecuaciones normales.
Resolviendo este sistema se obtienen los coeficientes buscados
m
!
∑ xi2
!
m
∑ yi −
!
Mm
!
∑ xiyi
m
∑ xi m
m
∑ xiyi
!
−
m
∑ xi
!
m
∑ yi
!
A
i=1 i=1 i=1 i=1 i=1 i=1 i=1
a0 = ! !2 , a1 = ! !2 .
m m m m
∑ xi2 − ∑ xi2 −
-F

m ∑ xi m ∑ xi
i=1 i=1 i=1 i=1

Para los datos del problema inicial se obtienen a0 = −0.360 y a1 = 1.538, y el gráfico de la
aproximación lineal puede verse en la Figura (3).
o
ic
ér
um
.N
A

Figura 3: Modelo lineal

Ahora consideremos el caso general donde se tienen los puntos {(xi , yi )} para i = 1, . . . , m y
se propone un modelo polinomial de grado n, con n < m − 1
n
Pn (x) = an xn + an−1 xn−1 + · · · + a1 x1 + a0 = ∑ a j xij .
j=0

3
 Al igual que en el caso de dos coeficientes, se deben determinar a0 , a1 , . . . , an que minimizan
m
E(a0 , . . . , an ) = ∑ [yi − Pn (xi )]2 .
i=1

Podemos reescribir esta función convenientemente

21
m
E(a0 , . . . , an ) = ∑ [yi − Pn(xi)]2
i=1
m m m

20
= ∑ y2i − 2 ∑ Pn(xi)yi + ∑ (Pn(xi))2
i=1 i=1 i=1
! !2
m m n m n
= ∑ y2i − 2 ∑ ∑ a j xij yi + ∑ ∑ a j xij

F
i=1 i=1 j=0 ! i=1 j=0 ! !
m m n m n n
= ∑ y2i − 2 ∑ yi ∑ a j xij +∑ ∑ a j xij ∑ ak xik

A
i=1 i=1 j=0 ! i=1 j=0 k=0 !
m n m n n m
=
i=1
M
∑ y2i − 2 ∑ a j ∑ yixij
j=0 i=1
+∑ ∑ a j ak ∑ xij+k
j=0 k=0 i=1
A
Para minimizar E se deben calcular las derivadas parciales ∂ E/∂ a j para j = 0, . . . , n e igualar
a cero:
-F
!
m n m
∂E j j+k
= −2 ∑ yi xi + 2 ∑ ak ∑ xi = 0, para j = 0, . . . , n.
∂aj i=1 k=0 i=1

Así tenemos las n + 1 ecuaciones normales en las n + 1 incógnitas a0 , . . . , an :

o

n m m
ic

j+k j
∑ ak ∑ xi = ∑ yi xi , para j = 0, . . . , n.
k=0 i=1 i=1
ér

Es decir, se obtiene el siguiente sistema lineal

 m m m m m
um

0 1 2



 a0∑ i x + a 1∑ i x + a 2 ∑ xi + ... + an ∑ xin = ∑ yixi0,

 i=1 i=1 i=1 i=1 i=1

 m m m m m
1 2
∑ xi3 + . . . + an ∑ xin+1 = ∑ yixi1,



 a0 ∑ i x + a 1 ∑ i x + a 2
i=1 i=1 i=1 i=1 i=1
.N

.. .. .. .. ..
. . . . .







 m m m m m
n n+1
+ a2 ∑ xin+2 + . . . + xi2n ∑ yixin.

 a0 ∑ xi + a1 ∑ xi an ∑ =


A

i=1 i=1 i=1 i=1 i=1

Este sistema lineal puede expresarse matricialmente, y es posible probar que si las xi son
todas distintas entonces las ecuaciones normales admiten una única solución.

Observación: a veces el modelo propuesto no es polinomial. Por ejemplo i) y = b eax o

ii) y = b xa , donde a y b son los coeficientes a determinar.
Lamentablemente, si se repitiera el procedimiento anterior cambiando el polinomio por los
modelos i) o ii) no se obtiene un sistema lineal, sino un sistema no lineal el cual no se puede
resolver por métodos sencillos.

4
 Un método más simple consiste en aplicar logaritmo natural en los modelos i) o ii):

ln y = lnb + a x o ln y = lnb + a ln x,

de esta manera es posible usar nuevamente un modelo lineal.

21
Ejemplo: supongamos que se conocen los siguientes datos

xi 1.00 1.25 1.50 1.75 2.00

20
yi 5.10 5.79 6.53 7.45 8.46

y se sabe que corresponden a un modelo de la forma y = beax . Aplicando logaritmo se obtiene

F
ln y = ln b + ax, un modelo lineal ỹ = b̃ + ax. Antes de aplicar las fórmulas de cuadrados

A
mínimos para el modelo lineal conviene reescribir la tabla anterior:

xi
ỹi M 1.00
1.629
1.25
1.756
1.50
1.876
1.75
2.008
2.00
2.135
A
Así se obtienen a = 0.5056 y b̃ = 1.122 y por lo tanto b = eb̃ = e1.122 = 3.071.
-F
o
ic
ér
um
.N
A

5
 Clase 11 - Aproximación de funciones (2)

Supongamos ahora que f ∈ C[a, b] y se desea determinar el mejor polinomio (en el sentido
de cuadrados mínimos) Pn (x) de grado ≤ n que minimice la siguiente medida del error entre
la función f y Pn en el intervalo [a, b]:

21
Z b Z b n
E = E(a0 , . . . , an ) = [ f (x) − Pn (x)]2 dx = [ f (x) − ∑ ak xk ]2 dx,
a a k=0

20
es decir, se deben determinar los coeficientes a0 , . . . , an que definen el polinomio Pn (x) de
manera que E sea mínima. Ver Figura 1.

F
A
M
A
-F
o
ic

Figura 1: Gráficos de f y Pn .
ér

Al igual que antes, una condición necesaria para encontrar un minimizador en a0 , . . . , an es

que ∂ E/∂ a j = 0 para todo j = 0, . . . , n. Antes de calcular estas derivadas, vamos a reescribir
convenientemente la expresión E:
um

Z b n
E = E(a0 , . . . , an ) = [ f (x) − ∑ ak xk ]2 dx
a k=0
Z b n Z b Z b n
.N

2
= [ f (x)] dx − 2 ∑ ak k
x f (x) dx + [ ∑ ak xk ]2 dx,
a k=0 a a k=0

luego para j = 0, . . . , n,
A

Z b n Z b
∂E j
= −2 x f (x) dx + 2 ∑ ak xk+ j dx = 0.
∂aj a k=0 a

Así, se obtienen las ecuaciones normales:

n Z b Z b
k+ j
∑ ak x dx = x j f (x) dx para j = 0, . . . , n.
k=0 a a

1
 Ejemplo: determinar el polinomio de aproximación de cuadrados mínimos de grado ≤ 2
para la función f (x) = sen π x en el intervalo [0, 1].
Las ecuaciones normales para el polinomio P2 (x) = a2 x2 + a1 x + a0 , están dadas por:

R1 R1 R1 2 R1
a0 0 1 dx + a1 0 x dx + a2 0 x dx = 0 sen πx dx

21
R1 R1 2 R1 3
R1
a0 0 x dx + a1 0 x dx + a2 0 x dx = 0 x sen πx dx
R1 2 R1 3 R1 4 R1 2
a0 0 x dx + a1 0 x dx + a2 0 x dx = 0 x sen πx dx

20
Calculando las integrales este sistema lineal puede escribirse matricialmente
    
1 1/2 1/3 a0 2/π

F
 1/2 1/3 1/4   a1  =  1/π 
1/3 1/4 1/5 a2 2
(π − 4)/π 3

A
de donde se obtiene que

a0 =
12π 2 − 120
π3 M
≈ −0.050465 a1 = −a2 =
720 − 60π 2
π3
≈ 4.12251
A
y el polinomio de grado ≤ 2 de mejor aproximación por cuadrados mínimos está dado por
P2 (x) = −4.12251 x2 + 4.12251 x − 0.050465. Ver Figura 2.
-F
o
ic
ér
um
.N

Figura 2: Gráficos de sen(x) y P2 .

Los coeficientes de la matriz de coeficientes pueden calcularse usando la siguiente fórmula

A

para el caso general

b j+k+1 − a j+k+1
Z b
a jk = x j+k dx = .
a j+k+1
Esta matriz es llamada matriz de Hilbert y es conocida por ser mal condicionada. El concep-
to de condicionamiento de una matriz se estudia en álgebra lineal numérica y está asociado
a la sensibilidad de la matriz frente a perturbaciones en los datos. Se dice que una matriz
es mal condicionada cuando pequeñas perturbaciones en los datos producen grandes pertur-
baciones en las soluciones. Esta es una característica de la matriz de coeficientes y no del
método numérico que se utilice para resolver el sistema lineal. Independientemente de esto,
sería deseable que la matriz sea lo más simple posible, por ejemplo una matriz diagonal.

2
 Para lograr esto, en lugar de proponer un polinomio de la forma Pn (x) = a0 + a1 x + · · · +
an xn , como una combinación lineal de los polinomios {x j }nj=0 , consideraremos otra forma
de expresar el mismo polinomio de cuadrados mínimos. A continuación veremos algunas
definiciones básicas y resultados que serán necesarios.
Definición 1. El conjunto de funciones {φ0 , . . . , φn } es linealmente independiente en el
intervalo [a, b], si siempre que

21
c0 φ0 (x) + c1 φ1 (x) + · · · + cn φn (x) = 0 para cualquier x ∈ [a, b],

se tiene que c0 = c1 = . . . cn = 0. En caso contrario se dice que ese conjunto de funciones es

20
linealmente dependiente.

F
Teorema 1. Si φ j (x) es un polinomio en x de grado igual a j para j = 0, . . . , n, entonces

A
{φ0 , . . . , φn } es un conjunto linealmente independiente para cualquier intervalo [a, b].

M
Demostración. Sean c0 , . . . , cn números reales tales que

P(x) = c0 φ0 (x) + c1 φ1 (x) + · · · + cn φn (x) =

n
∑ c j φ j (x) = 0,
A
j=0
-F

para cualquier x ∈ [a, b].

Como P(x) se anula en todo el intervalo [a, b], los coeficientes de todas las potencias de x
son iguales a cero. Como cn φn (x) es el único término que incluye xn entonces el coeficiente
n−1
o

cn = 0 y por lo tanto, ∑ c j φ j (x).

j=0
ic

Repitiendo esta misma idea, se tiene que el único término que incluye a xn−1 es cn−1 φn−1 (x)
y de aquí se concluye que cn−1 = 0. De igual forma se obtiene que cn−2 = · · · = c1 = c0 = 0,
ér

y, en consecuencia, {φ0 , . . . , φn } es un conjunto linealmente independiente.

um

El siguiente resultado para conjuntos de polinomios es análogo a otro que se utiliza en aplica-
ciones de álgebra lineal. No veremos la demostración, pero, en cambio, haremos un ejemplo.
.N

Teorema 2. Si {φ0 , . . . , φn } es un conjunto de polinomios linealmente independiente para

cualquier intervalo [a, b] en el espacio de polinomios de grado ≤ n, entonces todo polinomio
de grado ≤ n puede escribirse, de manera única, como combinación lineal de {φ0 , . . . , φn }.
A

Ejemplo: Sean φ0 (x) = 2, φ1 (x) = x − 3 y φ2 (x) = x2 + 2x + 7. Por el Teorema 1, {φ0 , φ1 , φ2 }

es un conjunto linealmente independiente en cualquier intervalo [a, b]. Sea Q(x) = a0 +a1 x +
a2 x2 un polinomio cuadrático arbitrario. Aplicaremos el teorema anterior mostrando que
existen coeficientes c0 , c1 y c2 tales que Q(x) = c0 φ0 (x) + c1 φ1 (x) + c2 φ2 (x).
1 3
Notar que 1 = φ0 (x), x = φ1 (x) + 3 = φ1 (x) + φ0 (x), y que
2 2
   
2 3 1 13
x = φ2 (x)−2x−7 = φ2 (x)−2 φ1 (x) + φ0 (x) −7 φ0 (x) = φ2 (x)−2φ1 (x)− φ0 (x).
2 2 2

3
 Luego

Q(x) = a0 + a1 x + a2 x2
= a0 12 φ0 (x) + a1 φ1 (x) + 23 φ0 (x) + a2 φ2 (x) − 2φ1 (x) − 13




2 φ0 (x)

= 12 a0 + 32 a1 − 13


2 a 2 φ0 (x) + (a1 − 2a2 )φ1 (x) + a2 φ2 (x),

21
y por lo tanto cualquier polinomio cuadrático puede escribirse como una combinación lineal
de {φ0 , φ1 , φ2 }.

20
F
A
M
A
-F
o
ic
ér
um
.N
A

4
 Clase 12 - Aproximación de funciones (3)

Mejor aproximación de funciones con pesos

Ahora se reformulará el problema de mejor aproximación de funciones por cuadrados míni-

mos con funciones de peso.

21
Definición 1. Una función (integrable) ω se llama función de peso en el intervalo I si
ω(x) ≥ 0 para todo x ∈ I, pero ω(x) 6= 0 para todo x en cualquier subintervalo de I, es decir

20
ω no puede ser constantemente cero en un subintervalo de I.

Las funciones de peso se utilizarán en la definición de la medida del error y permiten dar
más o menos importancia a las aproximaciones en ciertas partes del intervalo. Por ejemplo,

F
la función de peso definida por

A
1
ω(x) = √ para x ∈ (−1, 1)
1 − x2
M
pone menos énfasis cerca del origen y por el contrario mucho más cerca de los extremos del
intervalo. Ver Figura (1).
A
-F
o
ic

Figura 1: Función de peso w(x) = √ 1

1−x2
ér

El problema reformulado: Dados {φ0 , φ1 , . . . , φn } un conjunto de funciones linealmente

independientes en [a, b], una función de peso ω definida en [a, b] y f una función continua
um

en [a, b], se desean determinar los coeficientes a0 , a1 , . . . , an de la combinación lineal que

definen
n
P(x) = ∑ ak φk (x),
.N

k=0

tal que minimizan la medida del error

Z b Z b n
A

E = E(a0 , . . . , an ) = ω(x)[ f (x) − P(x)]2 dx = ω(x)[ f (x) − ∑ ak φk (x)]2 dx.

a a k=0

Notar que esto generaliza lo que vimos en la clase anterior tomando ω(x) ≡ 1 y φk (x) = xk
para k = 0, . . . , n.
Nuevamente, una condición necesaria para encontrar un minimizador en (a0 , . . . , an ) es que
∂ E/∂ a j = 0 para todo j = 0, . . . , n, esto es,
Z b n
∂E
=2 ω(x)[ f (x) − ∑ ak φk (x)] φ j (x) dx = 0, para j = 0, . . . , n.
∂aj a k=0

1
 Así, obtenemos las ecuaciones normales para el caso general
n Z b Z b
∑ ak ω(x)φk (x)φ j (x) dx = ω(x) f (x)φ j (x) dx, para j = 0, . . . , n.
k=0 a a

Si fuera posible elegir funciones {φ0 , . . . , φn } tal que

21
(
Z b 0 si j 6= k
ω(x)φk (x)φ j (x) dx = (1)
a αj > 0 si j=k

20
las ecuaciones normales se podrían reducir a
Z b Z b
aj ω(x)(φ j (x))2 dx = ω(x) f (x)φ j (x) dx, para j = 0, . . . , n.

F
a a

Ahora, usando (1) se tiene

A
Z b
a jα j = ω(x) f (x)φ j (x) dx, para j = 0, . . . , n,

y por lo tanto
a

M
A
Z b
1
aj = ω(x) f (x)φ j (x) dx, para j = 0, . . . , n.
αj a
-F

Definición 2. El conjunto {φ0 , . . . , φn } es un conjunto ortogonal de funciones en el intervalo

[a, b], con respecto a la función de peso ω, si
(
Z b 0 si j 6= k
o

ω(x)φk (x)φ j (x) dx =

a α j si j = k.
ic

Si α j = 1 para todo j = 0, . . . , n se dice que el conjunto es ortonormal.

ér

Lema 1. Si {φ0 , φ1 , . . . , φn } es un conjunto ortogonal de funciones en el intervalo I con

respecto a una función de peso ω definida en I entonces son linealmente independientes.
um

Demostración. Supongamos que

n
∑ c j φ j (x) = 0, para todo x ∈ I.
.N

j=0

Luego
A

Z b Z b n n Z b
0= 0 φk (x) ω(x) dx = ∑ c j φ j (x) φk (x) ω(x) dx = ∑ c j φ j (x) φk (x) ω(x) dx = ck αk ,
a a j=0 j=0 a

como αk > 0, entonces ck = 0 para todo k = 0, . . . , n, y por lo tanto las funciones son lineal-
mente independientes.
Lema 2. Sea el conjunto de funciones polinomiales {φ0 , φ1 , . . . , φn } es un conjunto ortogonal
en el intervalo [a, b] con respecto a una función de peso ω, con grado de φk igual a k y Qk (x)
es un polinomio de grado k menor estricto que n entonces
Z b
ω(x)φn Qk (x) dx = 0.
a

2
 Demostración. Como Qk (x) tiene grado k se sabe que existen coeficientes c0 , . . . , ck tales
que
k
Qk (x) = ∑ c j φ j (x).
j=0

Luego
Z b k Z b

21
ω(x)φn (x)Qk (x) dx = ∑ cj ω(x)φn (x)φ j (x) dx = 0,
a j=0 a

pues φn es ortogonal a φ j para cada j = 0, . . . , k.

20
En resumen, con lo anterior se ha probado el siguiente teorema

Teorema 1. Si {φ0 , φ1 , . . . , φn } es un conjunto ortogonal de funciones en el intervalo [a, b]

F
con respecto a una función de peso ω definida en [a, b], entonces la aproximación por cua-
drados mínimos a una función continua f respecto al peso ω está dada por

A
n
P(x) = ∑ ak φk (x),
donde para cada k = 0, . . . , n,
M k=0
A
Rb Z b
a ω(x) f (x)φk (x) dx
1
-F

ak = R b = ω(x) f (x)φk (x) dx.

2
a ω(x)(φk (x)) dx
αk a

El resultado siguiente da una relación de recurrencia que permite generar un conjunto de

o

funciones ortogonales.
ic

Teorema 2. El conjunto de funciones polinomiales {φ0 , φ1 , . . . , φn } que se define a conti-

nuación es un conjunto ortogonal en el intervalo [a, b] con respecto a una función de peso
ér

ω
φ0 (x) = 1, φ1 (x) = x − B1 para cada x ∈ [a, b],
um

donde Rb
xω(x)(φ0 (x))2 dx
B1 = Rab ,
2
a ω(x)(φ0 (x)) dx
y para k ≥ 2
.N

φk (x) = (x − Bk ) φk−1 (x) −Ck φk−2 (x) para cada x ∈ [a, b],
A

donde Rb Rb
a xω(x)(φk−1 (x))2 dx a xω(x)φk−1 (x)φk−2 (x) dx
Bk = R b y Ck = Rb .
a ω(x)(φk−1 (x))2 dx a ω(x)(φk−2 (x))2 dx

Demostración. La prueba se hará por inducción en k.

Si k = 1, entonces,
Z b Z b Z b
ω(x)φ1 (x)φ0 (x) dx = xω(x)φ0 (x) dx − B1 ω(x)φ0 (x) dx = 0,
a a a

pues φ0 (x) = (φ0 (x))2 y por la definición de B1 .

3
 Ahora supongamos, por hipótesis inductiva, que {φ0 , φ1 , . . . , φk−1 } son ortogonales, y vea-
mos que φk es ortogonal a todas las funciones anteriores.
Z b Z b
ω(x)φk (x)φk−1 (x) dx = ω(x) ((x − Bk ) φk−1 (x) −Ck φk−2 (x)) φk−1 (x) dx
a a
Z b Z b
2
= xω(x)(φk−1 (x)) dx − Bk ω(x)(φk−1 (x))2 dx
a a

21
Z b
−Ck ω(x)φk−2 (x)φk−1 (x) dx = 0,
a

20
pues los dos primeros términos suman cero por la definición de Bk y el último término es
cero por la hipótesis inductiva.
Además,

F
Z b Z b
ω(x)φk (x)φk−2 (x) dx = ω(x) ((x − Bk ) φk−1 (x) −Ck φk−2 (x)) φk−2 (x) dx

A
a a
Z b Z b
= xω(x)φk−1 (x)φk−2 (x) dx − Bk ω(x)φk−1 (x)φk−2 (x) dx
a

−CkM Z b

a
ω(x)(φk−2 (x))2 dx = 0,
a
A
pues el primero y el tercer término suman cero por la definición de Ck y el segundo término
-F

es cero por la hipótesis inductiva.

Por último para 0 ≤ i ≤ k − 3, reemplazando k por i + 1 en φk (x) = (x − Bk ) φk−1 (x) −
Ck φk−2 (x), se obtiene que φi+1 (x) = (x − Bi+1 ) φi (x) −Ci+1 φi−1 (x) y por lo tanto
o

xφi (x) = φi+1 (x) + Bi+1 φi (x) +Ci+1 φi−1 (x). (2)
ic

Luego
ér

Z b Z b
ω(x)φk (x)φi (x) dx = ω(x) ((x − Bk ) φk−1 (x) −Ck φk−2 (x)) φi (x) dx
a a
um

Z b Z b
= xω(x)φk−1 (x)φi (x) dx − Bk ω(x)φk−1 (x)φi (x) dx
a a
Z b
−Ck ω(x)φk−2 (x)φi (x) dx.
a
.N

Los dos últimos términos son iguales a cero por la hipótesis inductiva. Usamos (2) para
analizar el primer término:
A

Z b Z b
xω(x)φk−1 (x)φi (x) dx = ω(x)xφi (x)φk−1 (x) dx
a a
Z b
= ω(x) (φi+1 (x) + Bi+1 φi (x) +Ci+1 φi−1 (x)) φk−1 (x) dx
a
Z b Z b
= ω(x)φi+1 (x)φk−1 (x) dx + Bi+1 ω(x)φi (x)φk−1 (x) dx
a a
Z b
= +Ci+1 ω(x)φi−1 (x)φk−1 (x) dx,
a

y estos tres términos son iguales a cero por la hipótesis inductiva si 0 ≤ i ≤ k − 3. Por lo tanto
queda demostrada la ortogonalidad en todos los casos.

4
 Ejemplos:

- Polinomios de Legendre: I = [−1, 1], ω(x) = 1 para todo x ∈ I,

1 3 6 3
φ0 (x) = 1, φ1 (x) = x, φ2 (x) = x2 − , φ3 (x) = x3 − x, φ4 (x) = x4 − x2 + , . . .

21
3 5 7 25

- Polinomios de Chebyshev: I = (−1, 1), ω(x) = √ 1 para todo x ∈ I,

1−x2

20
φk+1 (x) = cos(k arc cos(x)), para todo x ∈ I.

- Polinomios de Laguerre: I = [0, +∞), ω(x) = e−x para todo x ∈ I,

F
ex d k −x k

A
φk+1 (x) = (e x ), para todo x ∈ I.
k! dxk

M 2
- Polinomios de Hermite: I = (−∞, +∞), ω(x) = e−x para todo x ∈ I,
A
2 d k x2
φk+1 (x) = (−1)k ex e , para todo x ∈ I.
dxk
-F
o
ic
ér
um
.N
A

5
 Clase 13 - Integración numérica
La integración numérica es una herramienta de gran utilidad en las siguientes aplicaciones:
- Se desea calcular la integral definida de f en el intervalo [a, b], pero f es muy complicada
o no tiene primitiva, como por ejemplo

21
Z 1
2
e−x dx
0

20
- no se tiene la función f en forma explícita sino sólo se conocen algunos valores funcionales,
por ejemplo una lista de pares ordenados (xi , yi ) para i = 0, . . . , n representados en un gráfico.
En ambos casos se desea estimar aproximadamente el valor de

F
Z b
f (x)dx.

A
a

Los métodos básicos que veremos son también conocidos como cuadratura numérica y
tienen la forma
M n
∑ ai f (xi).
i=0
A
Los métodos de cuadratura numérica se basan en interpolación numérica. Consideremos el
-F

conjunto de nodos distintos {x0 , . . . , xn } en el intervalo [a, b]. Sea Pn el polinomio, que inter-
pola a f en esos puntos, en la forma de Lagrange y en el término del error correspondiente
n
f n+1 (ξx ) n
Pn (x) = ∑ f (xi )Li (x), en (x) = (x − xi ),
(n + 1)! ∏
o

i=0 i=0
ic

para algún ξx ∈ (x0 , xn ). Es decir f (x) = Pn (x) + en (x). Luego, integrando

Z b Z b
ér

Rb
a f (x) dx = Pn (x) dx + en (x) dx
a a
Z b n
(ξx ) n
Z b (n+1)
f
um

=
a i=0
∑ f (xi)Li(x) dx + a (n + 1)! ∏(x − xi) dx
i=0
n Z b n
1 (n+1)
= ∑ ai f (xi ) + f (ξx ) ∏(x − xi ) dx
i=0 (n + 1)! a i=0
.N

Z b
para algún ξx ∈ (x0 , xn ) y donde ai = Li (x) dx para i = 0, . . . , n.
a
Por lo tanto las fórmulas de cuadratura numérica están dadas por
A

Z b n
f (x)dx ≈ ∑ ai f (xi ),
a i=0

con un error de integración numérica dado por

Z b n
1
En ( f ) = f (n+1) (ξx ) ∏(x − xi ) dx.
(n + 1)! a i=0

Cuando los puntos de interpolación estén igualmente espaciados, estas reglas o fórmulas se
llaman cuadratura de Newton–Cotes. Inicialmente veremos las reglas simples y luego las
reglas compuestas.

1
 Reglas simples

Las variantes que veremos a continuación dependen de la cantidad de puntos de interpolación

que se utilicen.

Regla del trapecio

21
Z b
Para integrar f (x) dx vamos a considerar dos puntos x0 = a, x1 = b y h = b − a = x1 − x0 .
a

20
Ver Figura 1. De esta manera se aproxima a la función f por el polinomio interpolante
x − x1 x − x0
P1 (x) = f (x0 ) + f (x1 ), (1)
x0 − x1 x1 − x0

F
y su correspondiente error es

A
f 00 (ξx )
e1 (x) = (x − x0 )(x − x1 ). (2)
2!

M
A
-F
o
ic

Figura 1: Regla del trapecio

Figura .5. Gráfica Regla del Trapecio Simple

Luego por (1),

ér

Z b Z b 
x − x1 x − x0
f (x) dx ≈ f (x0 ) + f (x1 ) dx
a a x0 − x1 x1 − x0
um

Z b Z b
x − x1 x − x0
= f (x0 ) dx + f (x1 ) dx.
a x0 − x1 a x1 − x0

(x1 − x0 ) b − a h
.N

Es muy fácil verificar que ambas integrales son iguales a = = , y por lo

2 2 2
tanto, Z b
(b − a) h
f (x) dx ≈ IT ( f ; a, b) = [ f (a) + f (b)] = [ f (a) + f (b)] .
2 2
A

Para poder expresar el error de la integral numérica de una manera más simple se utiliza
el siguiente teorema de Análisis Matemático, cuya demostración no se incluirá pero puede
consultarse en libros de Cálculo.
Teorema 1. Supongamos que f ∈ C[a, b], que g es una función integrable en [a, b] y que g
no cambia de signo en [a, b]. Entonces existe c ∈ (a, b) tal que
Z b Z b
f (x)g(x) dx = f (c) g(x) dx.
a a
Rb 1 Rb
En particular, si g(x) ≡ 1, entonces a f (x) dx = f (c)(b − a), esto es, f (c) = b−a a f (x) dx.

2
 Como g(x) = (x − x0 )(x − x1 ) = (x − a)(x − b) no cambia de signo en [a, b] y aplicando el
teorema anterior se sabe que existe ξ independiente de x tal que
Z b Z b
f 00 (ξx )(x − a)(x − b) dx = f 00 (ξ ) (x − a)(x − b) dx
a a
Z b
= f 00 (ξ ) (x2 − (a + b)x + ab) dx

21
a
b
x3 (a + b) 2

00
= f (ξ ) − x + abx . (3)

20
3 2 a

Luego,
 3 b
x (a + b) 2 b3 ab2 b3 a3 a3 a2 b

F
− x + abx = − − + ab2 − + + − a2 b
3 2 a 3 2 2 3 2 2

A
1 3
= (2b − 3ab2 − 3b3 + 6ab2 − 2a3 + 3a3 + 3a2 b − 6a2 b)
6

=M 1
6
(−b3 + 3ab2 + a3 − 3a2 b)
A
(a − b)3 (b − a)3 h3
=− =− .
-F

= (4)
6 6 6
Entonces, por (2), (3) y (4), el error de la fórmula del trapecio está dado por:
Z b
1 b 00
Z
ET = E1 (x) = e1 dx = f (ξx )(x − a)(x − b) dx
o

2! a
a 3 
h3 (b − a)3 00
ic

1 00 h
= f (ξ ) − = − f 00 (ξ ) = − f (ξ ),
2! 6 12 12
ér

para algún ξ ∈ (a, b).

Resumiendo, la regla del trapecio simple para integración numérica en el intervalo [a, b] está
um

dada por
Z b
(b − a) h
f (x) dx ≈ IT ( f ; a, b) = [ f (a) + f (b)] = [ f (a) + f (b)] ,
a 2 2
.N

y su correspondiente error es
(b − a)3 00 h3
ET = − f (ξ ) = − f 00 (ξ ),
12 12
A

para algún ξ ∈ (a, b).

Observación: como el término del error tiene f 00 , se dice que la regla del trapecio será
exacta cuando se aplica a una función tal que f 00 (x) ≡ 0, esto es, cualquier polinomio de
grado menor o igual que 1. Se puede verificar esta conclusión estudiando el comportamiento
de la regla del trapecio cuando se integran polinomios.
Z b
- Si f (x) ≡ 1, por un lado 1 dx = b − a,
a
(b − a)
y por otro, IT ( f ; a, b) = (1 + 1) = b − a. Luego, la regla del trapecio integra
2
exactamente cualquier polinomio de grado 0.

3
 b2 − a2
Z b
- Si f (x) ≡ x, por un lado x dx = ,
a 2
(b − a) b2 − a2
y por otro, IT ( f ; a, b) = (a + b) = . Luego, la regla del trapecio integra
2 2
exactamente cualquier polinomio de grado 1.

b3 − a3
Z b
- Si f (x) ≡ x2 , entonces x2 dx = ,

21
a 3
(b − a) 2 b3 + a2 b − ab2 − a3
y en cambio, IT ( f ; a, b) = (a +b2 ) = . Por lo tanto la regla
2 2

20
del trapecio no integra exactamente a polinomios de grado 2.

Regla de Simpson

F
Z b

A
En este caso, para integrar f (x) dx vamos a considerar tres puntos de interpolación x0 = a,
a
a+b b−a
x1 =
2
y x2 = b. Si llamamos h =
2
M entonces, x1 = x0 + h y x2 = a + 2h = b. Ver
Figura 2. Con tres puntos (n = 2) se construye un polinomio interpolante de grado 2, de
manera análoga a la aplicada en la regla del trapecio.
A
-F
o
ic
ér

Figura 2: Regla de Simpson

um

La deducción completa de la regla de Simpson es un ejercicio del práctico. El cálculo del

error es un poco más complicado y no será demostrado por falta de tiempo. En resumen, la
regla de Simpson está dada por:
.N

Z b  
h (b − a)/2 a+b
f (x) dx ≈ IS ( f ; a, b) = [ f (a) + 4 f (a + h) + f (a + 2h)] = f (a) + 4 f ( ) + f (b) ,
a 3 3 2

y su correspondiente error es
A

h5 (4) ((b − a)/2)5 (4)

ES = − f (ξ ) = − f (ξ ),
90 90

para algún ξ ∈ (a, b).

Observación: como el término del error tiene una derivada de orden 4, la regla de Simpson
integrará exactamente cuando f (4) ≡ 0, esto es, para polinomios de grado menor o igual a 3.

Por último veremos 2 reglas más sencillas que utilizan un único punto de interpolación.

4
 Regla del rectángulo

Esta regla utiliza sólo uno de los extremos de integración para construir el polinomio inter-
polante, por ejemplo se considera sólo x0 = a y por lo tanto el polinomio interpolante será
una constante (n = 0). Ver Figura 3.

21
20
F
A
Figura
Figura .2. 3: ReglaRegla
Gráfica del rectángulo
del Rectángulo

En este caso, la regla del rectángulo está dada por:

Z b
Mf (x) dx ≈ IR ( f ; a, b) = f (a)(b − a),
A
a

y su correspondiente error es
-F

(b − a)2 0
ER = f (ξ ),
2
para algún ξ ∈ (a, b).
Observación: también es posible tomar x0 = b, haciendo los cambios correspondientes en
o

la fórmula de IR , Además, como el término del error tiene una derivada de orden 1, la regla
ic

del rectángulo integrará exactamente a polinomios de grado 0.

ér

Regla del punto medio

um

Esta regla también utiliza sólo un punto pero en este caso es el punto medio del intervalo
x0 = (a + b)/2. Así, el polinomio interpolante será una constante (n = 0). Ver Figura 4.
.N
A

Figura 4: Regla del punto medio

Figura .3. Gráfica Regla del Punto Medio Simple

En este caso, la regla del punto medio está dada por:

Z b  
a+b
f (x) dx ≈ IPM ( f ; a, b) = f (b − a),
a 2

5
 y su correspondiente error es
(b − a)3 00
EPM = f (ξ ),
24
para algún ξ ∈ (a, b).
Observación: como el término del error tiene una derivada de orden 2, la regla del punto
medio integrará exactamente a polinomios de grado 1.

21
La siguiente definición es de fundamental importancia en el estudio de las reglas de integra-

20
ción numérica.

Definición 1. La precisión o grado de exactitud de una fórmula o regla de cuadratura es

el mayor entero no negativo n tal que la fórmula de integración es exacta para xk , para todo

F
k = 0, . . . n.

A
Así los métodos considerados en esta clase tienen la siguiente precisión.

M
Regla de cuadratura
Rectángulo
Precisión
0
A
Punto medio 1
Trapecio 1
-F

Simpson 3

En la siguiente tabla se resumen las reglas simples de integración numérica para estimar
o

Z b
f (x) dx:
ic

a
ér

Regla Puntos Fórmula Error Precisión

um

(b−a)2 0
Rectángulo 1 f (a)(b − a) 2 f (ξ ) 0
(b−a)3 00
Punto medio 1 f ( a+b
2 )(b − a) 24 f (ξ ) 1
.N

(b−a) 3
Trapecio 2 2 [ f (a) + f (b)] − (b−a) 00
12 f (ξ ) 1
(b−a)/2  5
f (a) + 4 f ( a+b − ((b−a)/2) f (4) (ξ )

Simpson 3 3 2 )+ f (b) 90 3
A

6
 Clase 14 - Integración numérica (2)
Las fórmulas de Newton-Cotes y las reglas simples de integración numérica en general no
son eficientes cuando se desea calcular una integral definida en un intervalo grande. Por un
lado, se requerirían más puntos de interpolación y los valores de los coeficientes son más
difíciles de calcular. Por otro lado, se sabe que los polinomios interpolantes suelen tener un

21
comportamiento oscilatorio cuando se utilizan muchos puntos de interpolación. Las reglas
compuestas que veremos a continuación se basan en particionar el intervalo de integración y
usar las reglas simples que ya vimos.

20
Reglas compuestas

F
Para entender mejor en que consisten las reglas compuestas se analizará un ejemplo aplican-
do la Regla de Simpson.

A
Z 4
Sea desea estimar el valor de ex dx. Esta integral es fácil de calcular con métodos analíti-
cos: Z 4
M 0

ex dx = [ex ]40 = e4 − e0 = e4 − 1 = 53.59815 . . .

A
0

Si aplicamos la Regla de Simpson

-F

Z 4
2
ex dx ≈ (e0 + 4e2 + e4 ) = 56.76958 . . .
0 3
esta estimación tendría un error aproximado de −3.17143 . . .
o

Ahora, si dividimos el intervalo en dos: [0, 4] = [0, 2] ∪ [2, 4], y aplicamos la Regla de Sim-
ic

pson en cada subintervalo se tiene

ér

Z 4 Z 2 Z 4
ex dx = ex dx + ex dx
0 0 2
1 0 1 + e2 ) + 1 (e2 + 4e3 + e4 )
≈ (e + 4e
um

3 3
1 0 1 2 3 4
≈ 3 (e + 4e + 2e + 4e + e ) = 53.8635 . . . ,

con una estimación del error de −0.26570 . . .

.N

Si dividimos otra vez el intervalo [0, 4] = [0, 1] ∪ [1, 2] ∪ [2, 3] ∪ [3, 4], y aplicamos la Regla
de Simpson en cada subintervalo se tiene
Z 4 Z 1 Z 2 Z 3 Z 4
A

x x x x
e dx = e dx + e dx + e dx + ex dx
0 0 1 2 3
1 0 1/2 + e) + 1 (e + 4e3/2 + e2 )
≈ 6 (e + 4e 6

+ 16 (e2 + 4e5/2 + e3 ) + 61 (e3 + 4e7/2 + e4 )

1 0 1/2 + 2e + 4e3/2 + 2e2 + 4e5/2 + 2e3 + 4e7/2 + e4 )
= 6 (e + 4e

= 53.61622 . . . ,

con una estimación del error de −0.01807 . . .

A continuación se generalizará esta idea para la Regla de Simpson.

1
 Regla compuesta de Simpson

Consideremos n par y subdividimos el intervalo [a, b] en n subintervalos iguales y aplique-

mos la regla de Simpson en cada subintervalo. Como n es par, se tiene una cantidad impar
de puntos igualmente espaciados x j = a + jh, para j = 0, . . . , n, con h = (b − a)/n.
Luego,

21
Z b n/2 Z x2 j
f (x) dx = ∑ f (x) dx
a j=1 x2 j−2

20
n/2 n   h5 (4)
h o
= ∑ f (x2 j−2 ) + 4 f (x2 j−1 ) + f (x2 j ) − f (ξ j ) ,
j=1 3 90

F
para algún ξ j ∈ (x2 j−2 , x2 j ), con j = 1, . . . , (n/2), y f ∈ C4 [a, b].
Notar que para j = 1, . . . , (n/2)−1, el término f (x2 j ) aparece en los subintervalos [x2 j−2 , x2 j ]

A
y [x2 j , x2 j+2 ], entonces:
Z b

a
f (x) dx =
hn
3 M
(n/2)−1 n/2 o h5 n/2
f (x0 ) + 2 ∑ f (x2 j ) + 4 ∑ f (x2 j−1 ) + f (xn ) −
j=1 j=1 90 ∑ f (4)(ξ j ),
j=1
A
para algún ξ j ∈ (x2 j−2 , x2 j ), con j = 1, . . . , (n/2), y f ∈ C4 [a, b].
-F

Consideremos el último término, correspondiente al error.

Como f (4) es continua en [a, b], entonces por el Teorema de valores extremos para funciones
continuas, se tiene que
o

mı́n f (4) (x) ≤ f (4) (ξ j ) ≤ máx f (4) (x), para j = 1, . . . , n/2

x∈[a,b] x∈[a,b]
ic

n/2
n n
mı́n f (4) (x) ≤ ∑ f (4)(ξ j ) ≤ máx f (4) (x)
ér

2 x∈[a,b] j=1 2 x∈[a,b]

n/2
2
mı́n f (4) (x) ≤ ∑ f (4)(ξ j ) ≤ máx f (4) (x).
um

x∈[a,b] n j=1 x∈[a,b]

Por el teorema del Valor Intermedio para funciones continuas, existe µ ∈ (a, b) tal que

2 n/2 (4)
.N

(4)
f (µ) = ∑ f (ξ j ),
n j=1

y por lo tanto,
A

n/2
n
∑ f (4)(ξ j ) = 2 f (4)(µ). (1)
j=1

Usando (1) y que h = (b − a)/n, el término del error en la regla compuesta de Simpson puede
ser reformulado independientemente de ξ j :

h5 n/2 (4) h5 (b − a) 4 (4)

E( f ) = − ∑ f (ξ j ) = − n f (4) (µ) = − h f (µ).
90 j=1 180 180

Los resultados desarrollados hasta aquí se resumen en el siguiente teorema:

2
 Teorema 1. Sean f ∈ C4 [a, b], n par, h = (b − a)/n y x j = a + jh, para j = 0, . . . , n. Entonces
existe µ ∈ (a, b) tal que la regla compuesta de Simpson para n subintervalos está dada por:

Z b (n/2)−1 n/2 o (b − a)
hn
f (x) dx = f (x0 ) + 2 ∑ f (x2 j ) + 4 ∑ f (x2 j−1 ) + f (xn ) − h4 f (4) (µ).
a 3 j=1 j=1 180

21
A continuación se presenta el pseudocódigo de la regla compuesta de Simpson.

20
Algoritmo de la regla compuesta de Simspon

F
A
Dados los siguientes datos de entrada: a: extremo inferior de integración, b: extremo superior
de integración y n un entero positivo par correspondiente al número de subintervalos en que
se particiona el [a, b].
input a, b, n M
A
h ← (b − a)/n
sx0 ← f (a) + f (b)
-F

sx1 ← 0 (suma de f (x2 j−1 ))

sx2 ← 0 (suma de f (x2 j ))
o

x←a
ic

for j = 1, 2, . . . , n − 1 do
x ← x+h
ér

if j es par then
um

sx2 ← sx2 + f (x)

else
sx1 ← sx1 + f (x)
.N

endif
endfor
sx ← (sx0 + 2 sx2 + 4 sx1) ∗ h/3
A

output sx
end

Regla compuesta del Trapecio

La deducción de esta regla se hace manera análoga a lo que se hizo con la regla compuesta
de Simpson. En este caso comenzaremos enunciando el siguiente teorema.

3
 Teorema 2. Sean f ∈ C2 [a, b], n un número entero positivo, h = (b − a)/n y x j = a + jh,
para j = 0, . . . , n. Entonces existe µ ∈ (a, b) tal que la regla compuesta del Trapecio para
n subintervalos está dada por:
Z b n−1 o (b − a)
hn
f (x) dx = f (a) + 2 ∑ f (x j ) + f (b) − h2 f 00 (µ).
a 2 j=1 12

21
Demostración. Se comienza particionando el intervalo [a, b] en n subintervalos y luego se
aplica la regla simple del Trapecio en cada subintervalo (Figura (1)):

20
Z b n Z xj
f (x) dx = ∑ f (x) dx
a j=1 x j−1
n n   h3 00
h

F
o
= ∑ f (x j−1 ) + f (x j ) − f (ξ j )
j=1 2 12

A
para algún ξ j ∈ (x j−1 , x j ), con j = 1, . . . , n, y f ∈ C2 [a, b].

M
Notar que los valores f (x j ) para j = 1, . . . , n − 1, aparecen dos veces en la última expresión,
entonces se puede escribir
A
Z b n−1 o h3 n
hn
f (x) dx = f (x0 ) + 2 ∑ f (x j ) + f (xn ) − ∑ f 00(ξ j ),
-F

a 2 j=1 12 j=1

para algún ξ j ∈ (x j−1 , x j )), con j = 1, . . . , n.

o
ic
ér
um

Figura 1: Regla compuesta del Trapecio

.N

Ahora, consideramos el término del error

n
h3
E( f ) = − ∑ f 00(ξ j ),
A

12 j=1

para algún ξ j ∈ (x j−1 , x j )), con j = 1, . . . , n. Como f 00 es continua en [a, b], entonces por el
Teorema de valores extremos para funciones continuas, se tiene que

mı́n f 00 (x) ≤ f 00 (ξ j ) ≤ máx f 00 (x), para j = 1, . . . , n/2

x∈[a,b] x∈[a,b]
n
n mı́n f 00 (x) ≤ ∑ f 00(ξ j ) ≤ n máx f 00 (x)
x∈[a,b] j=1 x∈[a,b]
n
1
mı́n f 00 (x) ≤ ∑ f 00(ξ j ) ≤ máx f 00 (x).
x∈[a,b] n j=1 x∈[a,b]

4
 Por el teorema del Valor Intermedio para funciones continuas, existe µ ∈ (a, b) tal que
n
1
f 00 (µ) = ∑ f 00(ξ j ),
n j=1

y por lo tanto,
n
∑ f 00(ξ j ) = n f 00(µ). (2)

21
j=1

Luego Usando (2) y que h = (b − a)/n, el término del error en la regla compuesta del Trape-

20
cio puede ser reformulado independientemente de ξ j :
n
h3 h3 (b − a) 2 (00)
E( f ) = −
12 ∑ f 00(ξ j ) = − 12 n f 00(µ) = − 12
h f (µ).
j=1

F
A
El pseudocódigo es muy fácil de deducir a partir la fórmula de la regla compuesta del Tra-
pecio.

M
Reglas compuestas del Punto Medio y del Rectángulo
A
Las deducciones de las reglas compuestas del Punto Medio y del Rectángulo son análogas
-F

a las dos anteriores por lo que sólo enunciamos los teoremas siguientes. Los respectivos
pseudocódigos se deducen fácilmente a partir de estos teoremas.
Teorema 3. Sean f ∈ C2 [a, b], n un número par, h = (b − a)/(n + 2) y x j = a + ( j + 1)h,
o

para j = −1, 0, . . . , n + 1. Entonces existe µ ∈ (a, b) tal que la regla compuesta del Punto
Medio para n + 2 subintervalos está dada por:
ic

Z b n/2
(b − a) 2 00
f (x) dx = 2h ∑ f (x2 j ) + h f (µ).
ér

a j=0 6

Ver Figura (2).

um
.N
A

Figura 2: Regla compuesta del Punto Medio

Teorema 4. Sean f ∈ C1 [a, b], n un número entero positivo, h = (b − a)/n y x j = a + jh, para
j = 0, . . . , n. Entonces existe µ ∈ (a, b) tal que la regla compuesta del Rectángulo para n
subintervalos está dada por:
Z b n−1
(b − a) 0
f (x) dx = h ∑ f (x j ) + h f (µ).
a j=0 2

5
 Ejemplo: si bien es posible calcular exactamente la integral
Z π
sen x dx = [− cos x]π0 = [cos x]0π = 2,
0

se desea determinar cuál es el entero n para estimar esta integral numéricamente con un error
menor que 0.00002, utilizando las diferentes reglas compuestas de integración.

21
- Regla compuesta del Rectángulo.
Usamos la expresión del error para despejar n y el hecho que h = (b − a)/n:

20
(b − a) 0
ER ( f ) = h f (µ),
2
luego,

F
(π − 0) π − 0 π2
|ER ( f )| = cos(µ) = < 0.00002,
2 n 2n

A
de aquí resulta que n > 246741.

M
- Regla compuesta del Punto Medio.
Usamos la expresión del error para despejar n y el hecho que h = (b − a)/(n + 2):
A
(b − a) 2 00
EPM ( f ) = h f (µ),
-F

6
luego,
2
π3

(π − 0) π −0
|EPM ( f )| = (− sen(µ)) = < 0.00002,
o

6 n+2 6(n + 2)2

ic

de aquí resulta que n > 507.

ér

- Regla compuesta del Trapecio.

Usamos la expresión del error para despejar n y el hecho que h = (b − a)/n:
um

(b − a) 2 00
ET ( f ) = − h f (µ),
12
luego,
.N

2
π3

(π − 0) π −0
|ET ( f )| = − (− sen(µ)) = < 0.00002,
12 n 12n3

de aquí resulta que n > 360.

A

- Regla compuesta de Simpson. Usamos la expresión del error para despejar n y el

hecho que h = (b − a)/n:
(b − a) 4 (4)
ES ( f ) = − h f (µ),
180
luego,
4
π5

(π − 0) π −0
|ET ( f )| = − sen(µ) = < 0.00002,
180 n 180n4
de aquí resulta que n > 18.

6
 En la siguiente tabla se resumen las reglas compuestas de integración numérica para estimar
Z b
f (x) dx:
a

21
Regla Fórmula Error
n−1
(b−a) 0
Rectángulo h ∑ f (x j ) 2 h f (µ)

20
j=0
n/2
(b−a) 2 00
Punto medio 2h ∑ f (x2 j ) 6 h f (µ)
j=0

F
n−1
hn o
Trapecio f (a) + 2 ∑ f (x j ) + f (b) − (b−a) 2 00
12 h f (µ)
2

A
j=1
(n/2)−1 n/2
hn o
Simpson f (x0 ) + 2 ∑ f (x2 j ) + 4 ∑ f (x2 j−1 ) + f (xn ) − (b−a) 4 (4) (µ)
180 h f
3
M j=1 j=1
A
-F
o
ic
ér
um
.N
A

7
 Clase 15 - Integración numérica (3)

Reglas gaussianas

Las fórmulas de cuadratura consideradas hasta ahora, simples y compuestas, se construyen

usando valores funcionales en puntos conocidos. Es decir, todas tienen la forma

21
Z b n
f (x) dx ≈ ∑ ai f (xi ), (1)
a i=0

20
y son exactas para polinomios de grado ≤ n. En esta fórmula los nodos x0 , x1 , . . . , xn son
conocidos a priori y los coeficientes a0 , a1 , . . . , an se determinan unívocamente de manera
que (1) sea una igualdad, para ciertos polinomios.

F
Por ejemplo, si usáramos la regla simple del Trapecio, cuyos nodos son los extremos del
intervalo:

A
M
A
-F

Figura 1: Regla del Trapecio (simple)

o

En cambio, si se pudieran elegir adecuadamente los dos nodos se podría mejorar la precisión:
ic
ér
um
.N

Figura 2: Regla con dos nodos

Recordemos que:
A

Regla n Número de puntos (n + 1) Precisión

Rectángulo 0 1 0
Punto medio 0 1 1
Trapecio 1 2 1
Simpson 2 3 3

Es decir que para (n + 1) puntos, la precisión de cada regla de cuadratura es (n + 1) o n.

Las reglas gaussianas permiten seleccionar convenientemente los nodos, además de los co-
eficientes, de manera óptima en el sentido que la regla de integración sea exacta para polino-
mios del grado más alto posible (precisión). Es decir que se deberán determinar los (n + 1)

1
 nodos x0 , . . . , xn y los (n + 1) coeficientes a0 , . . . , an de manera que la regla de cuadratura con
funciones de peso
Z b n
w(x) f (x) dx ≈ ∑ ai f (xi ), (2)
a i=0
sea exacta para polinomios de grado ≤ 2n + 1.
Antes de estudiar el caso general, veamos el siguiente ejemplo.

21
Ejemplo: tomando la función de peso w(x) = 1, determinar los nodos x0 y x1 y los coefi-
cientes a0 y a1 tal que la regla de cuadratura

20
Z 1
f (x) dx ≈ a0 f (x0 ) + a1 f (x1 ), (3)
−1

F
sea exacta para polinomios p de grado ≤ 2 · 1 + 1 = 3.

A
- Si grado(p) = 0, basta considerar p(x) ≡ 1, entonces

M 2=
Z 1

−1
1 dx = a0 + a1 . (4)
A
- Si grado(p) = 1, basta considerar p(x) = x, entonces
-F

Z 1
0= x dx = a0 x0 + a1 x1 . (5)
−1

- Si grado(p) = 2, basta considerar p(x) = x2 , entonces

o
ic

Z 1
2
= x2 dx = a0 x02 + a1 x12 . (6)
3 −1
ér

- Si grado(p) = 3, basta considerar p(x) = x3 , entonces

um

Z 1
0= x3 dx = a0 x03 + a1 x13 . (7)
−1
.N

Si a0 = a1 = 0, por (4) se tiene un absurdo, por lo tanto no pueden ser simultáneamente cero.
Supongamos ahora uno de los dos coeficientes es cero y el otro no, por ejemplo que a0 = 0
y a1 6= 0. Luego por (4) se deduce que a1 = 2, y por (5) se tiene que x1 = 0, y por (6) se
A

llega a un absurdo. Análogamente si suponemos que a0 6= 0 y a1 = 0. Por lo tanto ambos

coeficientes deben ser distintos de cero.
De la ecuación (5), si x0 = 0, y como a1 6= 0, resulta que x1 = 0 y por (6) se llega a un
absurdo. Por lo tanto x0 no puede ser 0 y por un razonamiento análogo x1 tampoco es 0.
Luego a0 , a1 , x0 y x1 son distintos de 0.
De (5) y (7) se tiene que
a0 x03 −a1 x13
= ,
a0 x0 −a1 x1
de donde se deduce que x02 = x12 y por lo tanto |x0 | = |x1 |. Ahora se tienen dos casos.

2
 Si x0 = x1 , por (4) y (5) se obtendría que x0 = 0 y por lo tanto es absurdo. Por lo tanto,
x0 = −x1 .
Usando√(6) se sabe√que (a0 + a1 )x02 = 2/3, y entonces x02 = 1/3 y de aquí se tiene que
3 3
x0 = − y x1 = .
3 3
Por último, usando (4) y (5) se obtiene que a0 = a1 = 1.

21
El teorema siguiente, debido a Gauss, indica como deben elegirse los (n + 1) nodos de ma-

20
nera de obtener una fórmula de cuadratura que sea exacta para polinomios de hasta grado
2n + 1. El resultado será enunciado y probado para una función de peso w general, pero
puede ser aplicado al caso simple donde w(x) ≡ 1.

F
Teorema 1. Sea w una función de peso positiva definida en [a, b] y q un polinomio no nulo de
grado exactamente (n + 1) que es ortogonal a todo polinomio p de grado ≤ n (con respecto

A
a w), es decir,
Z b

M a
q(x)p(x)w(x) dx = 0.

Si x0 , x1 , . . . , xn son las (n + 1) raíces de q, entonces la fórmula

(8)
A
Z b n
f (x)w(x) dx ≈ ∑ ai f (xi ) (9)
-F

a i=0
Z b n x−xj
con ai = w(x) ∏ dx, será exacta para todo polinomio f de grado ≤ 2n + 1.
a j=0 xi − x j
o

j6=i
ic

Demostración. Veamos primero que la fórmula de cuadratura es exacta para polinomios de

grado ≤ n.
ér

Sea f un polinomio de grado ≤ n, entonces f es el único polinomio que interpola los n + 1

n n x−x
j
um

nodos x0 , x1 , . . . , xn , y por lo tanto f se puede reescribir como f (x) = ∑ f (xi ) ∏ .

i=1 j=0 xi − x j
j6=i

Luego
Z b n Z b n x−xj n
.N

f (x)w(x) dx = ∑ f (xi ) w(x) ∏ dx = ∑ ai f (xi ).

a i=1 a j=0 xi − x j i=0
j6=i

Por lo tanto la regla de cuadratura es exacta para polinomios de grado ≤ n.

A

Ahora supongamos que f es un polinomio tal que n + 1 ≤ grado( f ) ≤ 2n + 1. Dividimos

f por el polinomio q y, por el algoritmo de la división de polinomios, se sabe que existen
polinomios p (cociente) y r (resto) tal que

f (x) = p(x)q(x) + r(x), (10)

con grado(p) ≤ n (pues grado( f ) ≤ 2n + 1) y grado(r) ≤ n o r(x) ≡ 0.

Luego,
Z b
q(x)p(x)w(x) dx = 0, (11)
a

3
 por hipótesis, pues el grado de q es n + 1.
Además, como xi es una raíz de q(x) para i = 0, . . . , n, se tiene que
f (xi ) = p(xi )q(xi ) + r(xi ) = r(xi ), para i = 0, . . . , n. (12)

Luego como grado(r) ≤ n, entonces

21
Z b n
r(x)w(x) dx = ∑ ai r(xi ). (13)
a i=0

20
Finalmente, por (10), (11), (12) y (13), se tiene que
Z b Z b Z b n n
f (x)w(x) dx = (p(x)q(x) + r(x))w(x) dx = r(x)w(x) dx = ∑ ai r(xi ) = ∑ ai f (xi ),
a a a i=0 i=0

F
que es lo que se quería probar.

A
El teorema anterior afirma que las raíces de los polinomios ortogonales son los nodos más
convenientes para tener más precisión en las reglas gaussianas. El siguiente resultado resume

M
algunas propiedades importantes de estos puntos. La demostración puede consultarse en los
libros de la Bibliografía.
A
Teorema 2. Sean w una función de peso en [a, b] y {φ0 , . . . , φn } un conjunto polinomios tales
que grado(φk ) = k y son w-ortogonales en [a, b] en el sentido que
-F

Z b
φi (x)φ j (x)w(x) dx = 0, si i 6= j.
a
Si x0 , . . . , xk−1 son las k raíces de φk entonces tales raíces son reales, simples y pertenecen
o

al intervalo abierto (a, b).

ic
ér

Observación: si la función de peso que se utiliza en el intervalo [−1, 1] es w(x) ≡ 1 se

obtienen los polinomios de Legendre:
um

1
P0 (x) ≡ 1, P1 (x) = x, P2 (x) = x2 − ,
3
3 6 3
P3 (x) = x3 − x, P4 (x) = x4 − x2 + .
.N

5 7 35
Para estos polinomios se conocen sus raíces así como los coeficientes ai de las fórmulas de
cuadratura:
A

grado del polinomio raíces (xi ) coeficientes (ai )

2 -0.5773502692 1.0000000000
0.5773502692 1.0000000000
3 -0.7745966692 0.555555556
0.0000000000 0.888888889
0.7745966692 0.555555556
4 -0.8611363116 0.3478548451
-0.3399810436 0.6521451549
0.3399810436 0.6521451549
0.8611363116 0.3478548451

4
 La siguiente muestra como es posible utilizar las reglas gaussianas en diferentes intervalos.
Observación: supongamos que se conoce una fórmula de integración numérica
Z 1 n
f (x) dx ≈ ∑ ai f (xi ),
−1 i=0

21
para valores determinados de los coeficientes ai , para i = 0, . . . , n, que dependen de los nodos
x0 , . . . , xn ∈ [−1, 1] y se desea obtener una fórmula de integración numérica para

20
Z b
f (x) dx.
a

Esto se puede hacer muy fácilmente haciendo el cambio de variables t : [a, b] −→ [−1, 1],

F
dado por t = αx + β , donde los coeficientes α y β se determinan fácilmente resolviendo el
sistema lineal

A

α a + β = −1
α b+β = 1

y por lo tanto t =
b−a
2
b−a
x−
a+b
b−a M
. Entonces x =
b−a
2
t+
a + b dx
2
y
dt
=
b−a
2
, es decir,
A
dx = dt. Luego,
2
-F
Z b Z 1   Z 1  
b−a a+b b−a b−a b−a a+b
f (x) dx = f t+ dt = f t+ dt.
a −1 2 2 2 2 −1 2 2

Ejemplo: se desea estimar el valor de la integral

o

Z 1.5
2
e−x dx
ic

usando las reglas gaussianas con la función de peso w(x) ≡ 1 y usando polinomios de grado
ér

2 y 3. Se sabe que esa integral es aproximadamente ≈ 0.1093643. Notar que si el grado es

m, entonces el índice de los nodos de la regla de cuadratura es n = m − 1, pues los subíndices
um

comienzan desde 0. Por lo tanto para grados m = 2 y 3 se tienen n = 1 y 2 y la precisión

(2n + 1) será 3 y 5, respectivamente.
Haciendo el cambio de variables, resulta
1 5
.N

y = 4 x − 5, x= y+ .
4 4
Luego,
Z 1.5 Z 1
−x2 1 2 /16
e−(y+5)
A

e dx = dy.
1 4 −1

Si el grado es m = 2, es decir (n = 1) con los puntos x0 y x1 , entonces

Z 1.5
2
e−x dx ≈ 0.1094003
1

Si el grado es m = 3, es decir (n = 2) con los puntos x0 , x1 y x2 , entonces

Z 1.5
2
e−x dx ≈ 0.1093642
1

5
 Clase 16 - Resolución de sistemas lineales

El problema

Se desea encontrar x1 , . . . , xn solución del siguiente sistema de ecuaciones lineales:

21


 a11 x1 + a12 x2 + a13 x3 + . . . + a1n xn = b1
 21 x1 + a22 x2 + a23 x3 + . . . + a2n xn = b2
a



a31 x1 + a32 x2 + a33 x3 + . . . + a3n xn = b3

20
 .. .. .. .. .. .. .. .. .. .. ..
. . . . . . . . . . .




an1 x1 + an2 x2 + an3 x3 + . . . + ann xn = bn


Es conveniente escribir este sistema en forma matricial:

F
    
a11 a12 a13 . . . a1n x1 b1

A
 a21 a22 a23 . . . a2n   x2   b2 
    
 a31 a32 a33 . . . a3n   x3 b3
= ,
  

 ..
 .
M
..
.
..
.
.. ..
..
an1 an2 an3 . . . ann

  ..
 .
xn
 
 
..

bn
.


A
es decir A x = b, con A ∈ Rn×n y b ∈ Rn . Así, el problema de esta unidad puede formularse
-F

de la siguiente manera:

Dada A ∈ Rn×n y b ∈ Rn hallar x∗ solución de A x = b.

o

Este es uno de los problemas más importantes en Análisis Numérico en general porque la
resolución de muchos problemas de diferentes disciplinas o aplicaciones termina en la reso-
ic

lución de sistemas lineales. Para lo cual es muy importante disponer de métodos numéricos
y computacionales que sean rápidos y eficientes.
ér

A continuación enunciamos un teorema básico de Álgebra lineal que será de utilidad en lo

que sigue.
um

Teorema 1. para toda matriz A ∈ Rn×n , las siguientes propiedades son equivalentes:

1. existe A−1 , la inversa de A. (A−1 A = AA−1 = I);

.N

2. la matriz A es no singular, es decir, si Ax = 0 entonces x = 0;

3. det(A) 6= 0;
A

4. las columnas de A forman una base de Rn ;

5. las filas de A forman una base de Rn ;

6. para cada b ∈ Rn existe un único x ∈ Rn solución de A x = b.

Si bien se conocen métodos matemáticos de Álgebra lineal para resolver ese sistema lineal,
nuestro objetivo será estudiar métodos y algoritmos numéricos para resolver este problema
de una forma sistemática y eficiente. Esta eficiencia estará asociada al costo computacional
de un algoritmo. Este costo se mide en términos de cantidad de operaciones y en espacio
de almacenamiento. Dado que los problemas que estudiaremos son simples y pequeños nos

1
 ocuparemos sólo de contabilizar la cantidad de operaciones que realiza cada algoritmo para
obtener la solución buscada.
Antes de estudiar el caso de una matriz general, consideraremos la resolución de sistemas
lineales de dos casos muy simples.

Sistemas diagonales

21
Definición 1. Una matriz A es diagonal si y sólo si ai j = 0 para i 6= j.

20
Sea A ∈ Rn×n una matriz diagonal
 
a11 0 ... 0
 0 a22 ... 0 

F
A= .
 
.. .. .. ..
 . . . . 

A
0 0 . . . ann
Luego A es no singular si y sólo sí aii 6= 0, i = 1, . . . , n si y sólo si det(A) = a11 a22 . . . ann 6= 0.

M
Entonces Ax = b si y sólo si aii xi = bi , para i = 1, . . . , n, si y sólo si xi =
bi
aii
, para i = 1, . . . , n.
A
Algoritmo 1: resolución de sistema lineal con matriz diagonal.
-F

input n, A, b
for i = 1, . . . , n do
xi ← bi /aii
o

output i, xi
ic

end for
ér

end
um

Costo computacional: n operaciones (productos) en punto flotante, y así se tiene O(n) flops.

Observación: notar que como la matriz A es diagonal y no singular, los únicos elementos no
.N

nulos están precisamente en la diagonal. Entonces, sería suficiente almacenar los elementos
diagonales de A en un vector d = diag(A), y así calcular x(i) = b(i)/d(i) para i = 1, . . . , n.
A

Sistemas triangulares

Definición 2. Una matriz A es triangular inferior (superior) si y sólo si ai j = 0 para i < j

(ai j = 0 para i > j).

A ∈ Rn×n es triangular inferior y no singular si y sólo si det(A) = a11 a22 . . . ann 6= 0, con
 
a11 0 . . . 0
 a21 a22 . . . 0 
A= . .
 
. .. . .
 . . . 0 
an1 an2 . . . ann

2
 Entonces resolver Ax = b es equivalente a


 x1 = b1 /a11

 1
x2 = (b2 − a21 x1 )
 

 a11 x1 = b1 

 a22
.. ..
 

 a21 x1 + a22 x2 = b2 




 .. ..



 . .
. . i−1
1
 

21
=⇒ xi = (bi − ∑ ai j x j )
 ai1 x1 + · · · + aii xi = bi  aii
  j=1

 .. .. 

.. ..
. .

 

. .
 

20
 
an1 x1 + · · · + ann xn = bn
 

 n−1

 1

 xn = (b n − ∑ an j x j )
ann


j=1

F
Algoritmo 2: resolución de sistema lineal con matriz triangular inferior.

A
input n, A, b
for i = 1, . . . , n do

xi ←
i−1
bi − ∑ ai j x j /aii
!
M (no sumar si i = 1)
A
j=1

output i, xi
-F

end for
end
Costo computacional: vamos a contar las sumas/restas y productos/ divisiones que se reali-
o

zan para cada i y luego sumamos:

ic

i sumas productos
ér

1 0 1
2 1 2
3 2 3
um

.. .. ..
. . .
n n-1 n
(n − 1)n n(n + 1)
Total
.N

2 2

(n − 1)n n(n + 1)
Así, la cantidad total de operaciones es: + = n2 , es decir O(n2 ) flops.
2 2
A

n
n(n + 1)
(Recordar que ∑i= 2
).
i=1

Análogamente, Si A ∈ Rn×n es triangular superior y no singular, es decir,

 
a11 a12 . . . a1n
 0 a22 . . . a2n 
A= . ..  ,
 
. .. . .
 . . . . 
0 0 . . . ann
con det(A) = a11 a22 . . . ann 6= 0.

3
 Entonces resolver Ax = b es equivalente a


 xn = bn /ann
1


xn−1 = (bn−1 − an−1,n xn )
 


 a11 x1 + a12 x2 + · · · + a1n xn = b1 


 an−1,n−1


 a22 x2 + · · · + a2n xn = b2 

 .. ..

21


 .. ..


 . .

. .
 n
=⇒ 1
aii xi + · · · + ain xn = bi xi = (bi − ∑ ai j x j )

 
 aii j=i+1
.. ..

20

 


 . .

 .. ..
. .

 

ann xn = bn
 
 n
1


= (b1 − ∑ a1 j x j )

 x1


a11 j=2

F
Algoritmo 3: resolución de sistema lineal con matriz triangular superior.

A
input n, A, b
for i = n, . . . , 1 do
n
! M
A
xi ← bi − ∑ ai j x j /aii (no sumar si i = n)
j=i+1
-F

output i, xi
end for
end
o

Costo computacional: de manera análoga al caso anterior se puede ver que la cantidad total
de operaciones es: n2 , es decir O(n2 ) flops.
ic
ér

Eliminación gaussiana
um

Dado un sistema lineal Ax = b, la idea consiste en transformarlo en otro sistema equivalente

Ux = y, que sea más fácil de resolver. Recordar que sistemas equivalentes tienen la misma
solución x. La matriz U del nuevo sistema será triangular superior, cuya resolución ya vimos
que es muy sencilla. Recordemos también que existen 3 tipos de operaciones elementales por
.N

filas que permiten transformar un sistema lineal ampliado en otro sistema lineal ampliado
equivalente:

- reemplazar una fila por un múltiplo escalar (no nulo) de ella misma ( fi ← λ fi , λ 6= 0).
A

- restar a una fila un múltiplo escalar (no nulo) de otra fila ( fi ← fi − λ f j , λ 6= 0).
- intercambiar dos filas ( fi ↔ f j ).

El método de eliminación gaussiana se basa fundamentalmente en las dos primeras operacio-

nes elementales por filas. Es decir, se realizan (n − 1) pasos para obtener la matriz triangular
U partiendo de A:
Paso 1 Paso 2 Paso 3 Paso n − 1
A −−−→ A(1) −−−→ A(2) −−−→ . . . −−−−−→ A(n−1) = U
y las mismas operaciones que se realizan en A se deben realizar en el vector b.

4
 Para fijar ideas, veamos con detalles cómo se hacen estos pasos en un ejemplo numérico.
Consideremos el sistema lineal:

    
6 −2 2 4 x1 12
 12 −8 6 10   x2   34 
  = 
 3 −13 9 3   x3   27 
−6 4 1 −18 −38

21
x4

Paso 1: para pasar de la matriz A a la matriz A(1) se deben realizar operaciones elementales

20
por filas para obtener ceros en las posiciones de la primera columna a partir de la segunda
fila hasta la última.
Para esto, denotemos fi a la i-ésima fila de A. Llamamos a f1 como la fila pivote y a a11 = 6

F
será llamado el elemento pivote. Entonces haremos las siguientes operaciones:

A
f2 ← f2 − 2 f1
f3 ← f3 − (1/2) f1

Mf4 ← f4 − (−1) f1
A
Así obtenemos,     
6 −2 2 4 x1 12
-F

 0 −4 2 2   x2   10 
  = 
 0 −12 8 1   x3   21 
0 2 3 −14 x4 −26
o

Paso 2: para pasar de A(1) a A(2) se deben realizar operaciones elementales por filas para
obtener ceros en las posiciones de la segunda columna a partir de la tercera fila hasta la
ic

(1)
última. La fila pivote será la nueva f2 y el elemento pivote será a22 = −4. Se realizan las
siguientes operaciones:
ér

f3 ← f3 − 3 f2
f4 ← f4 − (−1/2) f2
um

Así obtenemos,    

6 −2 2 4 x1 12
 0 −4 2 2   x2   10 
=
.N

  
 0 0 2 −5   x3   −9 
0 0 4 −13 x4 −21

Paso 3: por último, para pasar de A(2) a A(3) = U se deben realizar operaciones elementales
A

por filas para modificar la última fila:

f4 ← f4 − 2 f3

Así obtenemos,    
6 −2 2 4 x1 12
 0 −4 2 2   x2   10 
  = 
 0 0 2 −5   x3   −9 
0 0 0 −3 x4 −3
cuya solución es: x = (1, −3, −2, 1).

5
 Observación 1: notar que, en general, al pasar de A a A(1) en el Paso 1, se obtiene una matriz
de la forma  (1) (1) (1)

a11 a12 . . . a1n
 
a11 a12 . . . a1n
(1) 
 a21 . . . a2n  0 . . . a2n 

 Paso 1 
..  −−−→ 

 .. .. . . .. .. . . .. 
 . . . .  
 . . .

. 
an1 an2 . . . ann 0 a
(1)
... a
(1)

21
n2 nn

Más específicamente, se calcula



20
 ai j si i = 1, j = 1, . . . , n 
 bi si i = 1



ai j = 0
(1) si i = 2, . . . , n, j = 1 y bi =
(1)
ai1
 ai1  bi − b1 si i = 2, . . . , n.
 ai j −

 a1 j si i = 2, . . . , n, j = 2, . . . , n a11

F
a11

Los coeficientes mi,1 = ai1 /a11 para i = 2, . . . , n se denominan multiplicadores.

A
Costo computacional: calcularemos por separado la cantidad de sumas/restas y de produc-
tos/divisiones.
Para el Paso 1: M
A
- Sumas/restas: (n − 1)2 (submatriz (n − 1) × (n − 1)) y (n − 1) (vector b).
-F

En total son: (n − 1)2 + (n − 1) = (n − 1)(n − 1 + 1) = n(n − 1) = n2 − n.

- Productos/divisiones: (n − 1) (multiplicadores), (n − 1)2 (submatriz (n − 1) × (n − 1))

y (n − 1) (vector b).
En total son: (n − 1)2 + 2(n − 1) = (n − 1)(n − 1 + 2) = (n − 1)(n + 1) = n2 − 1.
o
ic

Observación 2: para el Paso 2, es análogo a lo anterior, pero se trabaja con una matriz de
un orden menor, es decir (n − 1) × (n − 1), y por lo tanto se requieren: (n − 1)2 − (n − 1)
ér

sumas/restas y (n − 1)2 − 1 productos/divisiones.

Observación 3: en general, para pasar de la matriz A(k−1) a la matriz A(k) en el Paso k, se
um

calcula
 (k−1)


 ai j si i = 1, . . . , k, j = 1, . . . , n

 0 si i = k + 1, . . . , n, j = k

.N

(k)
ai j =
(k−1)

 (k−1) aik (k−1)
a − a si i = k + 1, . . . , n, j = k + 1, . . . , n


 ij
 (k−1) k j
akk
A

y
 (k−1)
b si i = 1, . . . , k
 i


(k) (k−1)
bi = (k−1) aik (k−1)
 bi

 − (k−1) bk si i = k + 1, . . . , n.
akk

Ejercicio: mostrar que el costo total de operaciones del método de eliminación gaussiana
para resolver un sistema Ax = b, con A ∈ Rn×n es O( 23 n3 ) flops.
n n
n(n + 1) n(n + 1)(2n + 1)
Ayuda: recordar que ∑ k = y ∑ k2 = .
k=1 2 k=1 6

6
 Algoritmo 4: eliminación gaussiana para obtener un sistema lineal con matriz triangular
superior.
input n, A, b
for k = 1, . . . , n − 1 do
for i = k + 1, . . . , n do

21
if (akk = 0) STOP!
m ← aik /akk

20
for j = k + 1, . . . , n do
ai j ← ai j − mak j

F
endfor ( j)
bi ← bi − mbk

A
endfor (i)
endfor (k)
output A, b
M
A
end
-F

Observación 1: notar que estamos almacenando la matriz triangular superior U en la misma

matriz A y el vector modificado y se sobreescribe en el mismo vector b. No es necesario
anular los elementos que están debajo de la diagonal, sólo se tiene en cuenta que se utilizarán
o

los elementos de la parte triangular superior.

ic

Observación 2: una vez obtenida la matriz triangular superior U y el vector y, se puede

aplicar el algoritmo 3 para resolver sistemas con matrices triangulares superiores.
ér

Observación 3: si para algún k se tiene que akk = 0, se podría intercambiar la fila k por otra
fila j, con j = k + 1, . . . , n tal que a jk 6= 0. Esto se llama estrategia de pivoteo y no se estudia
um

en este curso.
Observación 4: notar que el ciclo ( j) del algoritmo, donde se van recorriendo las columnas
de la matriz, comienza desde k +1 y no desde k. Esto se debe a que, en cada fila, el coeficiente
.N

que multiplica a la fila pivote (multiplicador) se define de manera de anular el coeficiente de

la columna k. Entonces no tiene sentido hacer las operaciones en la columna k pues ya se
sabe que el resultado va a dar cero.
A

Por último enunciamos un resultado que da condiciones suficientes para poder aplicar elimi-
nación gaussiana para resolver un sistema lineal.

Teorema 2. Sean Ak = A(1 : k, 1 : k) es la submatriz de A formada por las primeras k filas y

k columnas de A para todo k = 1, . . . , n. Si det(Ak ) 6= 0 para k = 1, . . . , n, entonces es posible
realizar el proceso completo de eliminación gaussiana y, por lo tanto, el sistema Ax = b tiene
única solución.

7
 Clase 17 - Resolución de sistemas lineales (2)

El problema

Dada A ∈ Rn×n y b ∈ Rn hallar x∗ solución de A x = b.

21
Factorización LU

20
Vimos que, usando eliminación gaussiana, esto se puede resolver con un costo computacional
de O( 32 n3 ) flops. Este costo computacional se debe principalmente a las operaciones que se
deben realizar para transformar la matriz A en una matriz triangular superior U, que es mucho

F
mayor que la cantidad de operaciones que se realizan para transformar el vector b en otro
vector y.

A
Ahora supongamos que se deben resolver varios sistema lineales Ax = b, con la misma matriz
A y diferentes vectores b. Sería poco eficiente repetir todas las operaciones elementales por

M
fila del proceso de eliminación gaussiana siendo que sólo el vector b es modificado y la
matriz A se mantiene igual.
A
La idea de la factorización LU consiste en factorizar la matriz como un producto de dos
matrices más simples: A = LU, donde L es una matriz triangular inferior con elementos
-F

diagonales iguales a uno (lii = 1, i = 1, . . . , n) y U triangular superior. Así:


Ly = b
Ax = b ⇐⇒ LUx = b ⇐⇒
Ux = y
o

Es decir que, conocidas las matrices L y U de la factorización de A, en lugar resolver Ax = b,

ic

se deben resolver dos sistemas triangulares para cada nuevo vector b.

Estas matrices L y U se obtienen a partir de A de una forma constructiva, como veremos a
ér

continuación. La idea es escribir A = LU, o equivalentemente LU = A, y tratar de despejar

los coeficientes de L y de U, usando inducción:
um

    
1 0 0 ... 0 u11 u12 u13 . . . u1n a11 a12 a13 . . . a1n
 l21 1 0 . . . 0   0 u22 u23 . . . u2n   a21 a22 a23 . . . a2n 
    
 l31 l32 1 . . . 0   0 0 u . . . u   a31 a32 a33 . . . a3n 
  33 3n = .
 .. .. .. .. ..   .. .. .. .. ..   .. .. .. .. .. 
.N

 . . . . .  . . . . .   . . . . . 
ln1 ln2 ln3 . . . 1 0 0 0 . . . unn an1 an2 an3 . . . ann

Vamos a determinar la primera fila de U y la primera columna de L.

A

Si multiplicamos la primera fila de L por la columna j de U para j = 1, . . . , n, obtenemos:

a1 j = 1 u1 j + 0 u2 j + 0 u3 j + · · · + 0 un j ,

por lo tanto,
u1 j = a1 j , para j = 1, . . . , n.

Si multiplicamos cada fila de L por la primera columna de U, obtenemos:

ai1 = li1 u11 , para i = 2, . . . , n,

1
 por lo tanto,
ai1
li1 = , para i = 2, . . . , n.
u11
Ahora, supongamos conocidas las (k − 1) filas de U y (k − 1) columnas de L y vamos a
determinar la fila k de U y la columna k de L.
Si multiplicamos la fila k de L por la columna j de U para j = k, . . . , n, obtenemos:

21
k−1
ak j = ∑ lkmum j + 1 uk j ,
m=1

20
por lo tanto,
k−1
uk j = ak j − ∑ lkmum j , para j = k, . . . , n.
m=1

F
Si multiplicamos cada fila de L, desde la k + 1 hasta la última, por la columna k de U,

A
obtenemos:
k−1
aik = ∑ limumk + lik ukk ,
por lo tanto,
1
M
k−1
m=1
A
lik = (aik − ∑ lim umk ) para i = k + 1, . . . , n.
ukk m=1
-F

A continuación veremos el algoritmo para realizar la factorización LU.

Algoritmo: Factorización LU.
o

input n, A = (ai j )
ic

for k = 1, . . . , n do
for j = k, . . . , n do
ér

k−1
uk j ← ak j − ∑ lkmum j
um

m=1

end for ( j)
for i = k + 1, . . . , n do (no ejecutar si k = n)
.N

k−1
1
lik ← (aik − ∑ lim umk )
ukk m=1

end for (i)

A

end for (k)

output L = (li j ),U = (ui j )
end
El siguiente resultado da una condición suficiente para poder realizar la factorización LU.
Teorema 1. Sea A ∈ Rn×n tal que las submatrices Ak = A(1 : k, 1 : k) para k = 1, . . . , n−1 son
no singulares. Entonces existen únicas matrices L,U ∈ Rn×n tal que L es triangular inferior
con 1 en la diagonal (lii = 1, i = 1, . . . , n) y U es triangular superior. Además, det(A) =
u11 u22 . . . unn .

2
 Demostración. La prueba se hará por inducción en la dimensión de la matriz.
Si n = 1, entonces si A = LU, basta tomar L = [1] y U = [a11 ]. También se cumple que
det(A) = u11 .
Supongamos que la factorización es válida para dimensión k − 1, es decir, dada Ak−1 ∈
R(k−1)×(k−1) , existen Lk−1 ,Uk−1 ∈ R(k−1)×(k−1) tal que Ak−1 = Lk−1Uk−1 . Veamos para k, o
sea, veremos que para una matriz Ak ∈ Rk×k existen Lk ,Uk ∈ Rk×k tal que Ak = LkUk . Para

21
determinar esto vamos a particionar convenientemente:
    
Ak−1 a Lk−1 0 Uk−1 e

20
= ,
bT c dT 1 0T ukk

donde a, b, d, e ∈ Rk−1 y c, ukk ∈ R.

F
Realizando los productos en bloques del lado derecho e igualando al lado izquierdo, obtene-
mos:

A
M Ak−1 = Lk−1Uk−1
Lk−1 e = a
(1)
(2)
A
d T Uk−1 = bT (3)
-F

T
d e + ukk = c (4)

De (2), y como Lk−1 es triangular inferior con unos en la diagonal, es claro que existe un
único e ∈ Rk−1 .
o

De (3), aplicando transpuesta a ambos lados, (d T Uk−1 )T = b, es decir, Uk−1

T d = b, y co-
ic

mo Uk−1 tiene inversa por hipótesis inductiva y por que las matrices Ak son no singulares,
entonces existe un único d ∈ Rk−1 .
ér

Luego, de (4), ukk = c − d T e está bien definido y está unívocamente determinado.

De esta manera se ha probado que existen Lk y Uk tales que Ak = LkUk .
um

Por último, como A = LU, entonces

det(A) = det(LU) = det(L) det(U) = 1 · (u11 u22 . . . unn ) = u11 u22 . . . unn .
.N
A

Observación 1: El costo de realizar la factorización LU es el mismo que para realizar la

eliminación gaussiana, es decir, O( 23 n3 ).
Observación 2: Si la matriz A no será usada posteriormente, una implementación eficiente
de la factorización LU, podría sobreescribir las matrices L y U sobre la misma A, teniendo en
cuenta que los coeficientes de la diagonal de L son todos iguales a 1. Más aún, se puede ver
que los elementos de la matriz L debajo de la diagonal son los multiplicadores que aparecían
en la eliminación gaussiana.
Observación 3: este es el método más eficiente para calcular el determinante de una matriz.
El costo computacional será aproximadamente el mismo que para realizar la factorización
LU.

3
 Recordemos la clásica regla para calcular el determinante para una matriz 2x2:
 
a b
det(A) = det = ad − bc
c d

Si usáramos esta regla se requieren 2 productos y una resta para calcular el determinante. La-
mentablemente, esta cantidad aumenta increíblemente cuando crece el orden n de la matriz.

21
Por ejemplo, si A fuera una matriz 20 × 20. Entonces se requieren del orden de 20! ≈ 2.4 1018
operaciones. Si se usara una computadora que realiza 109 operaciones en punto flotante por
segundo, entonces se necesitarían 76 años para calcular ese determinante.

20
Observación 4: nunca se debe calcular la inversa de una matriz. Siempre es mucho más
eficiente reformular el problema para resolver sistemas lineales. Sin embargo, si el único ob-
jetivo del problema es obtener la matriz de una inversa de una matriz, la forma más eficiente

F
consiste en realizar una factorización LU de la matriz A y luego resolver n sistemas lineales
de la forma LUx = ei , donde ei es el i-ésimo vector de la base canónica que tiene un 1 en la

A
posición i, y 0 en las otras componentes, para i = 1, . . . , n.

M
A continuación veremos algunas definiciones y resultados que serán de utilidad para la clase
siguiente.
A
Definición 1. Una norma vectorial en Rn es una función que asigna a cada x ∈ Rn un
-F

número real no negativo denotado por kxk y llamado norma de x, tal que se satisfacen las
siguientes tres propiedades para todo x, y ∈ Rn y para todo α ∈ R:

1. kxk > 0 si x 6= 0, y k0k = 0;

o

2. kαxk = |α|kxk;
ic

3. kx + yk ≤ kxk + kyk.
ér

Ejemplos:
um

!1/2
n
2
1. norma euclideana (o norma 2): kxk2 = ∑ |xi| ;
i=1
.N

n
2. norma 1: kxk1 = ∑ |xi |;
i=1
A

!1/p
n
3. para p ≥ 1, norma p: kxk p = ∑ |xi| p ;
i=1

4. norma infinito: kxk∞ = máx |xi |;

1≤i≤n

Definición 2. Dados dos vectores x, y ∈ Rn y una norma vectorial k · k se define la distancia

entre x e y por:
d(x, y) = kx − yk.

El concepto de normas puede extenderse al conjunto de matrices.

4
 Definición 3. Una norma matricial en Rn×n es una función que asigna a cada A ∈ Rn×n
un número real no negativo denotado por kAk y llamado norma de A, tal que se satisfacen
las siguientes cuatro propiedades para todo A, B ∈ Rn×n y para todo α ∈ R:

1. kAk > 0 si A 6= 0, y k0k = 0;

2. kαAk = |α|kAk;

21
3. kA + Bk ≤ kAk + kBk;

4. kABk ≤ kAkkBk.

20
Ejemplo: norma de Frobenius
!1/2

F
n n
kAkF = ∑ ∑ |ai j |2 .

A
i=1 j=1

Definición 4. Dada una norma vectorial en Rn y una matriz A ∈ Rn×n se define la norma
matricial inducida por
M kAk = sup
kAxk
x6=0 kxk
.
A
Proposición 1. Sean A ∈ Rn×n y k · k una norma vectorial en Rn que induce una norma
-F

matricial, entonces:

1. kAxk ≤ kAkkxk para todo x ∈ Rn ;

o

2. existe x̄ con kx̄k = 1 tal que kAx̄k = kAk.

ic

Demostración. - Probemos la primera parte.

Si x = 0, entonces kAxk = 0 = kAkkxk.
ér

kAxk kAxk
Si x 6= 0, entonces kAk = sup ≥ para todo x, y por lo tanto kAxk ≤ kAkkxk
x6=0 kxk kxk
um

para todo x ∈ Rn .

- Ahora probemos la segunda parte.

 
kAxk Ax x
.N

kAk = sup = sup = sup A = sup kAyk.

x6=0 kxk x6=0 kxk x6=0 kxk y:kyk=1

Como S = {x ∈ Rn |kxk = 1} es cerrado y acotado (compacto), la función continua

A

x → kAxk alcanza sus valores extremos, en particular existe un x̄ tal que kx̄k = 1 y

kAk = sup kAyk = máx kAyk = kAx̄k.

y:kyk=1 y:kyk=1

Ejemplos:

n
1. kAk1 = máx ∑ |ai j |;
1≤ j≤n i=1

5
 n
2. kAk∞ = máx ∑ |ai j |;
1≤i≤n j=1

q
3. kAk2 = ρ(AT A),
donde ρ(B) es el radio espectral de B y se define por máx{|λ | : λ es autovalor de B}.

21
20
F
A
M
A
-F
o
ic
ér
um
.N
A

6
 Clase 18 - Resolución de sistemas lineales (3)

El problema

Dada A ∈ Rn×n y b ∈ Rn hallar x∗ solución de A x = b.

Hasta ahora vimos dos métodos numéricos para resolver este problema: eliminación gaus-

21
siana y factorización LU. Ambos métodos son considerados métodos directos y se sabe que
luego de un número finito de pasos se obtiene una solución, salvo errores de redondeo. Por

20
otro lado, existe otra familia de métodos para resolver este problema llamados iterativos o
indirectos.

F
Métodos iterativos

A
Los métodos iterativos para sistemas lineales generan una sucesión de vectores {x(k) }, a par-
tir de un vector inicial x(0) , que convergen a la solución de Ax = b, bajo adecuadas hipótesis.

M
La convergencia significa que esta sucesión se detiene cuando se alcanza una precisión deter-
minada luego de un cierto número de iteraciones. El éxito computacional de estos métodos
requiere que procedimiento sea simple y de bajo costo computacional. En general estos mé-
A
todos son adecuados para matrices grandes y que tienen muchos coeficientes iguales a cero
-F

(matrices ralas).
Estudiaremos dos métodos iterativos muy conocidos: método de Jacobi y método de Gauss-
Seidel. Presentaremos brevemente estos métodos y algunos resultados simples de conver-
gencia. Para fijar ideas veremos inicialmente un ejemplo muy sencillo.
o

Ejemplo: consideremos el sistema lineal

ic

     
7 −6 x1 3 7x1 − 6x2 = 3
= ←→ , (1)
ér

−8 9 x2 −4 −8x1 + 9x2 = −4

cuya solución es x∗ = (0.2, −0.2666 . . . ). A partir del sistema lineal (1), podemos despejar
um

las variables:
x1 = 76 x2 + 37
(
(2)
x2 = 98 x1 − 49
.N

(k−1) (k−1)
Dada una aproximación x(k−1) = (x1 , x2 ), y de (2) se puede generar el siguiente mé-
todo iterativo: 
 x1(k) = 6 (k−1)
+ 37
7 x2
A

,
(k) 8 (k−1)
x2 = 9 x1 − 49


el cual es conocido como método de Jacobi.

También se puede generar otro método iterativo de la siguiente manera:

 x1(k) = 6 x2(k−1) + 3
7 7
,
 (k) 8 (k) 4
x2 = 9 x1 − 9

el cual es conocido como método de Gauss-Seidel.

1
 La tabla siguiente muestra algunas iteraciones de ambos métodos, comenzando en ambos
(0) (0)
casos con x(0) = (x1 , x2 ) = (0, 0). Recordemos que la solución es: x∗ = (0.2, −0.2666 . . . ).

It. M. de Jacobi M. de Gauss-Seidel

(k) (k) (k) (k)
k x1 x2 x1 x2
0 0 0 0 0

21
10 0.14865 -0.19820 0.21978 -0.24909
20 0.18682 -0.24909 0.20130 -0.26531
30 0.19662 -0.26215 0.20009 -0.26659

20
40 0.19913 -0.26551 0.20001 -0.26666
50 0.19978 -0.26637 0.20000 -0.26667

F
Los métodos de Jacobi y Gauss-Seidel son conocidos también como métodos de separa-
ción (splitting) y parten de una misma idea básica, que consiste en escribir la matriz como

A
A = M − N, donde M es una matriz no singular. Así tenemos que

M Ax = b
A
(M − N)x = b
Mx = Nx + b (3)
-F

x = M −1 (Nx + b) (4)
−1 −1
x = (M N)x + M b. (5)
o

(0) (0)
Dado x(0) = (x1 , x2 ), la última ecuación sugiere definir un método iterativo de la forma:
ic

x(k+1) = (M −1 N)x(k) + M −1 b, para k ≥ 0. (6)

ér

A continuación veremos algunos resultados de convergencia de estos métodos iterativos.

Teorema 1. Sea b ∈ Rn y A = M − N ∈ Rn×n donde A y M son matrices no singulares. Si
um

kM −1 Nk < 1 para alguna norma matricial inducida entonces la sucesión generada por la
ecuación (6) converge a la solución de Ax = b para cualquier vector inicial x(0) .

Demostración. Restando la ecuación (6) de (5), donde en (5) ponemos la solución x∗ se

.N

obtiene:
x(k+1) − x∗ = (M −1 N)(x(k) − x∗ ), para k ≥ 0. (7)
Ahora, utilizando alguna norma matricial inducida, se obtiene:
A

kx(k+1) − x∗ k ≤ k(M −1 N)kkx(k) − x∗ k, para k ≥ 0. (8)

Luego, reptiendo este último paso se tiene:
kx(k+1) − x∗ k ≤ k(M −1 N)kk+1 kx(0) − x∗ k, para k ≥ 0. (9)
Usando que k(M −1 N)k < 1, se tiene que

lı́m kx(k) − x∗ k = 0,
k→∞

para cualquier vector inicial x(0) .

2
 La principal diferencial entre el método de Jacobi y el de Gauss-Seidel consiste en cómo es
elegida la matriz M. Por simplicidad vamos a descomponer la matriz A de la siguiente forma:
A = L + D +U, donde L es la parte triangular inferior de A (sin la diagonal), D es la diagonal
de A y U es la parte triangular superior de A (sin la diagonal). Es importante no confundir
con las matrices L y U de la factorización LU.

21
Método de Jacobi

En el método de Jacobi se toma M = D, y por lo tanto,

20
N = M − A = D − (L + D +U) = −(L +U),

esto es,

F
   
a11 0 ... 0 0 a12 . . . a1n

A
 0 a22 ... 0   a21 0 . . . a2n 
M= y N = − ..  .
   
.. .. .. ..  .. .. ..
 . . . .   . . . . 
0 0
M
. . . ann

Si miramos la i-ésima componente en (3) tenemos que

an1 an2 ... 0
A
[Dx(k+1) ]i = [b − (L +U)x(k) ]i
-F


(k)
x1
(k+1)  . 
aii xi  .. 
= bi − [ai1 . . . ai,i−1 0 ai,i+1 . . . ain ]  
(k)
xn
o

i−1 n
(k+1) (k) (k)
aii xi = bi − ∑ ai j x j − ∑ ai j x j
ic

j=1 j=i+1
i−1 n
(k+1) 1 (k) (k)
ér

xi = (bi − ∑ ai j x j − ∑ ai j x j ),
aii j=1 j=i+1
um

para i = 1, . . . , n.
Ahora veremos un pseudocódigo del método de Jacobi. Notar que en el algoritmo usaremos
x para representar el vector inicial y sobreeescribiremos las sucesivas aproximaciones en el
mismo vector. La variable MaxIt representa el número máximo de iteraciones permitidas y
.N

Xtol la tolerancia para la norma de la diferencia de dos aproximaciones vectoriales sucesivas.

Algoritmo: método de Jacobi.
input n, A, b, x, MaxIt, Xtol
A

for k = 1, . . . , MaxIt do
for i = 1, . . . , n do
i−1 n
1
ui ← (bi − ∑ ai j x j − ∑ ai j x j )
aii j=1 j=i+1

end for (i)

if ku − xk < Xtol STOP, output: u es la solución, k iteración.
for i = 1, . . . , n do

3
 xi ← ui
end for (i)
end for (k)
output x
end

21
Antes de enunciar un resultado de convergencia, daremos una definición.
Definición 1. Una matriz A, de orden n × n, es diagonalmente dominante si

20
n
|aii | > ∑ |ai j |, para i = 1, . . . , n.
j=1
j6=i

F
Teorema 2. Si A es diagonalmente dominante, entonces la sucesión generada por el método

A
de Jacobi converge a la solución de Ax = b, para cualquier vector inicial x(0) ∈ Rn .

M
Demostración. En el método de Jacobi, la matriz M es la diagonal de A y debe ser inversible
para que el método esté bien definido, por lo que aii 6= 0, i = 1, . . . , n. La matriz de iteración
está dada por
A
  
1/a11 0 ... 0 0 a12 . . . a1n
-F

 0 1/a22 . . . 0   a21 0 . . . a2n 

 
M −1 N = −  .

.. .. . . .
..   .. . . . . . . .. 
 ..
 
. . 
0 0 . . . 1/ann an1 an2 . . . 0
  (10)
0 a12 /a11 . . . a1n /a11
o

 a21 /a22 0 . . . a2n /a22 

ic

= −
 
.. .. .. .. 
 . . . . 
an1 /ann an2 /ann . . . 0
ér

Luego,
n |ai j |
um

kM −1 Nk∞ = máx ∑ |aii| < 1,

1≤i≤nj=1
pues A es diagonalmente dominante. Finalmente, la convergencia es una consecuencia direc-
ta del teorema anterior.
.N

Método de Gauss-Seidel
A

En el método de Gauss-Seidel se toma M = L + D, y por lo tanto,

N = M − A = L + D − (L + D +U) = −U,
esto es,
   
a11 0 ... 0 0 a12 . . . a1n
 a21 a22 ... 0   0 0 . . . a2n 
M= y N = − ..  .
   
.. .. .. ..  .. .. ..
 . . . .   . . . . 
a1n a2n . . . ann 0 0 ... 0

4
 Luego por (3) tenemos que
(L + D)x(k+1) = b −Ux(k)
Dx(k+1) = b − Lx(k+1) −Ux(k) .
Así, analizando la i-ésima componente, obtenemos
[Dx(k+1) ]i = [b − Lx(k+1) −Ux(k) ]i

21
   
(k+1) (k)
x1 x1
(k+1)
aii xi

= bi − [ai1 . . . ai,i−1 0 . . . 0]  ..  − [0 . . . 0 ai,i+1 . . . ain ]  .. 
  
.  . 

20
 
(k+1) (k)
xn xn
i−1 n
(k+1) (k+1) (k)
aii xi = bi − ∑ ai j x j − ∑ ai j x j
j=1 j=i+1

F
i−1 n
(k+1) 1 (k+1) (k)
xi = (bi − ∑ ai j x j − ∑ ai j x j ),

A
aii j=1 j=i+1

para i = 1, . . . , n.

M
Ahora veremos un pseudocódigo del método de Gauss-Seidel. Notar que en el algoritmo usa-
remos x para representar el vector inicial y sobreeescribiremos las sucesivas aproximaciones
A
en el mismo vector. La variable MaxIt representa el número máximo de iteraciones permi-
tidas y Xtol la tolerancia para la norma de la diferencia de dos aproximaciones vectoriales
-F

sucesivas.
Algoritmo: método de Gauss-Seidel.
input n, A, b, x, MaxIt, Xtol
o

for k = 1, . . . , MaxIt do
ic

for i = 1, . . . , n do
ér

ui ← 0
end for (i)
um

for i = 1, . . . , n do
i−1 n
1
ui ← (bi − ∑ ai j u j − ∑ ai j x j )
aii j=1 j=i+1
.N

end for (i)

if ku − xk < Xtol STOP, output: u es la solución, k iteración.
A

for i = 1, . . . , n do
xi ← ui
end for (i)
end for (k)
output x
end
El resultado de convergencia que enunciamos para este método es similar al del método de
Jacobi, aunque su demostración no es tan directa.

5
 Teorema 3. Si A es diagonalmente dominante, entonces la sucesión generada por el método
de Gauss-Seidel converge a la solución de Ax = b, para cualquier vector inicial x(0) ∈ Rn .

Observación: existen resultados que muestran que si ambos métodos convergen, el méto-
do de Gauss-Seidel lo hace más rápido que el método de Jacobi. Esto es razonable, pues
al calcular la i-ésima componente, en el método de Gauss-Seidel, se utiliza información re-
cientemente calculada en la misma iteración. Por otro lado, si ambos métodos divergen, el

21
método de Gauss-Seidel lo hace más rápido que el método de Jacobi. Por esta razón el méto-
do era más atractivo computacionalmente que el de Jacobi. Sin embargo, el método de Jacobi

20
ha vuelta ha cobrar popularidad en las últimas décadas, pues es paralelizable mientras que el
método de Gauss-Seidel no lo es.
Por último vamos a enunciar un resultado auxiliar y un teorema que, bajo ciertas hipótesis,
asegura la convergencia de los métodos como Jacobi o Gauss-Seidel.

F
Teorema 4. Para cada matriz A ∈ Rn×n , se cumple que ρ(A) es igual al ı́nf{kAk} sobre

A
todas las normas matriciales inducidas.

M
Teorema 5. Una condición necesaria y suficiente para que la sucesión generada por el
método iterativo x(k+1) = (M −1 N)x(k) + M −1 b, para k ≥ 0, converja a la única solución de
Ax = b para todo x(0) inicial, es que ρ(M −1 N) < 1.
A
-F
o
ic
ér
um
.N
A

6
 Clase 19 - Programación lineal

Introducción

La programación lineal (PL) es uno de los mecanismos más naturales para tratar una gran
cantidad de problemas del mundo real de un modo sencillo. Es una subárea de Optimización

21
donde, como es de esperar, todas las funciones involucradas para formular el problema son
lineales. Aunque las funciones lineales son de las más simples, existen muchos problemas
de economía, producción, planning, logística, redes, scheduling, transporte y otras áreas que

20
pueden formularse como un problema de programación lineal. Además, los aspectos mate-
máticos y computaciones de PL son muy interesantes. Por un lado, la matemática de PL es
sencilla, poderosa y elegante y utiliza enfoques geómetricos y de Álgebra lineal, que están
conectados entre sí. Por otro lado, los aspectos computacionales son muy importantes por

F
las potenciales aplicaciones y variantes de los algoritmos para resolver el problema. El mé-

A
todo más conocido, es el Simplex, y fue propuesto por Dantzig en 1947. En estas clases, por
una cuestión de tiempo, presentaremos algunas nociones básicas de PL y el método Simplex
en una versión resumida. Para entender en que consiste un problema de PL consideraremos
inicialmente un ejemplo.
M
A
Ejemplo: un agricultor debe comprar fertilizantes (abono) para sus campos. El ingeniero
-F

agrónomo le dijo que cada kilogramo de fertilizante le alcanza para 10 m2 de su campo, y

debido a las características propias de esas tierras, el fertilizante debe contener (al menos):
3 g de fósforo (P), 1.5 g de nitrógeno (N) y 4 g de potasio (K) por cada 10 m2 . En el mercado
existen 2 tipos de fertilizantes: T1 y T2. El fertilizante T1 contiene 3 g de P, 1 g de N y 8 g
o

de K y cuesta $ 10 por kilogramo. En cambio, el fertilizante T2 contiene 2 g de P, 3 g de N

y 2 g de K y cuesta $ 8 por kilogramo. El agricultor desea saber cuántos kilogramos de cada
ic

fertilizante debe comprar, por cada 10 m2 de campo, de modo de minimizar el costo total
cubriendo los requerimientos de su suelo.
ér

Vamos a resumir toda esta información en la siguiente tabla y a continuación definiremos el

problema:
um

Tipo 1 (x1 ) Tipo 2 (x2 ) Necesidades mínimas

P (fósforo) 3 2 3
N (nitrógeno) 1 3 1.5
.N

K (potasio) 8 2 4
Costo 10 8

Incógnitas:
A

x : cantidad de fertilizante de tipo 1 (T1), en kg.

y : cantidad de fertilizante de tipo 2 (T2), en kg.
Restricciones (requerimientos):

3x + 2y ≥ 3
x + 3y ≥ 1.5
8x + 2y ≥ 4
x ≥ 0, y ≥ 0 (no negatividad)
Las últimas restricciones de negatividad son razonables pues la cantidad de fertilizantes no
puede ser negativa.

1
 Función objetivo (costo): f (x, y) = 10x + 8y.
Objetivo: Minimizar f (x, y) = 10x + 8y.
Así, este problema puede ser formulado en la siguiente forma:

Minimizar f (x, y) = 10x + 8y

sujeto a 3x + 2y ≥ 3

21
x + 3y ≥ 1.5
8x + 2y ≥ 4
x ≥ 0, y ≥ 0

20
o de una forma más compacta:

Minimizar cT x

F
sujeto a Ax ≥ b
x≥0

A
donde c = (10, 8), x = (x1 , x2 ) y

3 2
M


A= 1 3  y

3


b =  1.5 
A
8 2 4

Este problema se puede representar gráficamente como se ve en la Figura 1. La región som-

-F

breada se llama región factible y establece el conjunto de posibles soluciones, es decir los
posibles valores de x e y que satisfacen las restricciones.
o
ic
ér
um
.N

Figura 1: Representación gráfica del ejemplo de los fertilizantes.

A

El objetivo será encontrar la solución óptima en la región factible, es decir (x∗ , y∗ ) en la

región factible que minimiza la función objetivo (costo) f .
En general, los problemas de PL están definidos por:

- un vector de variables x ∈ Rn que son no negativas, es decir, xi ≥ 0 para i = 1, . . . , n.

- una función objetivo lineal f , que deberá ser minimizada o maximizada. Como esta
función debe ser lineal será de la forma f (x) = cT x = c1 x1 + · · · + cn xn .

- un conjunto de restricciones que deben satisfacer las variables. Estas restricciones es-
tarán dadas por ecuaciones (igualdades) o inecuaciones (desigualdades) lineales.

2
 - una región (o conjunto) factible definida por la intersección de todas las restricciones
lineales que deben satisfacer las variables.

Por lo tanto, los problemas de PL se escriben de la siguiente forma:

Minimizar f (x)
sujeto a Cx = d

21
(1)
Rx ≥ s
x≥0

20
donde f es la función objetivo a minimizar (o maximizar); el sistema lineal Cx = d co-
rresponde al conjunto de ecuaciones lineales de igualdad que deben satisfacer las variables;
Rx ≥ s corresponde al conjunto de restricciones lineales de desigualdad que deben satisfacer

F
las variables (podrían ser con ≤ en vez de ≥); y por último, las restricciones de no negativi-
dad x ≥ 0. Observar que en el ejemplo anterior no había restricciones de igualdad.

A
Se denomina forma estándar cuando el problema de PL es formulado como

M Minimizar
sujeto a
cT x
Ax = b
x≥0
(2)
A
-F

Veamos que cualquier problema de PL en formato general puede ser llevado a la forma
estándar:

- Conversión de maximización en minimización: si el problema original fuera:

o

maximizar z = 4x1 − 3x2 + 6x3 = cT x,

ic

multiplicando sólo la función objetivo por (−1) se convierte en

ér

minimizar ẑ = −4x1 + 3x2 − 6x3 = −cT x.

um

Una vez que el problema de minimización es resuelto, el valor de la función objetivo

se debe multiplicar por (−1), de manera que z∗ = −ẑ∗ , aunque las variables óptimas
son los mismas.
.N

- Si alguna componente bi < 0 para algún i, se debe multiplicar por (−1) a toda la
restricción, cambiando el signo de la desigualdad si la hubiera.

- Si la cota inferior de una variable no fuera cero, por ejemplo x1 ≥ 5, se puede definir
A

x̂1 = x1 − 5 ≥ 0.

- Si una variable tuviera una cota superior, por ejemplo x1 ≤ 7, se la tratará como cual-
quiera de las otras restricciones.

- Si una variable fuera irrestricta (libre), por ejemplo la variable x2 , entonces debe re-
emplazarse por: x2 = x̂2 − x20 , donde x̂2 , x20 ≥ 0.

- Una restricción de desigualdad con ≤ se convierte en una restricción de igualdad agre-

gando una variable de holgura (slack), por ejemplo: 2x1 + 7x2 − 3x3 ≤ 10 es equiva-
lente a 2x1 + 7x2 − 3x3 + s1 = 10 con s1 ≥ 0.

3
 - Una restricción de desigualdad con ≥ se convierte en una restricción de igualdad agre-
gando una variable de exceso, por ejemplo: 6x1 − 2x2 + 4x3 ≥ 15 es equivalente a
6x1 − 2x2 + 4x3 − s1 = 15 con s1 ≥ 0.

- Una restricción de igualdad puede convertirse en 2 restricciones de desigualdad, por

ejemplo: 2x1 + 3x2 = 1 es equivalente a 2x1 + 3x2 ≤ 1 y 2x1 + 3x2 ≥ 1.

21
Entonces, el problema del ejemplo puede ser escrito en la forma estándar:

Minimizar z = 10x1 + 8x2

20
sujeto a 3x1 + 2x2 − s1 = 3
x1 + 3x2 − s2 = 1.5
8x1 + 2x2 − s3 = 4
x1 , x2 , s1 , s2 , s3 ≥ 0

F
A
Ahora veremos algunas definiciones y resultados sobre la geometría del problema de PL.
Definición 1. Dados a = (a1 , a2 , . . . , an ) ∈ Rn , a 6= 0, y b ∈ R, se definen:

M
1. {x ∈ Rn |a1 x1 + a2 x2 + · · · + an xn = b} es llamado hiperplano (afín si b 6= 0) del
espacio n-dimensional ;
A
2. {x ∈ Rn |a1 x1 + a2 x2 + · · · + an xn ≤ b} es llamado semiespacio (cerrado) del espacio
-F

n-dimensional.

Observación: un hiperplano divide al espacio en dos semiespacios.

o

Ejemplos:
ic

1. x1 + x2 = 2 es un hiperplano en R2 ;
ér

2. x1 + x2 + x3 = 1 y x3 = 2 son 2 hiperplanos afines en R3 .

Definición 2. Un conjunto S de Rn es convexo si para todo par de puntos distintos x, y ∈ S,
um

el segmento que los une también está contenido en S, es decir, (αx + (1 − α)y) ∈ S para todo
α ∈ [0, 1]. (Ver Figura 2).
.N
A

Figura 2: Conjuntos convexo y no convexo.

Los siguientes dos lemas son muy fáciles de demostrar y sus demostraciones se dejan como
ejercicio.
Lema 1. la intersección finita de conjuntos convexos es un conjunto convexo.
Lema 2. todo semispacio cerrado es un convexo.

4
 Definición 3. La intersección de una cantidad finita de semiespacios cerrados de Rn se
denomina región poliedral cerrada de Rn .

Observación 1: por los dos lemas anteriores, toda región poliedral cerrada es un conjunto
convexo.
Observación 2: el conjunto de restricciones de un problema de programación lineal dado
por Ω = {x ∈ Rn |Ax = b, Rx ≤ s} es una región poliedral cerrada. Usualmente, el conjunto

21
factible de un problema de PL suele llamarse politopo, siempre que sea no vacío. Además, si
esta región es acotada suele llamarse siemplemente poliedro. Que sea acotada, significa que

20
si (x1 , . . . , xn ) ∈ Ω entonces existen constantes positivas ki tales |xi | ≤ ki para i = 1, . . . , n.

Para caracterizar las regiones poliedrales como Ω, estamos interesados en los vértices, los

F
cuales se determinan por la intersección de las ecuaciones que definen los semiespacios de
Ω.

A
En el Ejemplo 1, cuya región factible Ω está dada por las restricciones

M
3x + 2y
x + 3y
8x + 2y
≥
≥
≥
3
1.5
4
(r1)
(r2)
(r3)
A
x ≥ 0 (r4)
y ≥ 0 (r5)
-F

y está graficada en la Figura 1, los vértices son:

1 6 6 3 3
P1 = (0, 2), P2 = ( , ), P3 = ( , ), P4 = ( , 0).
5 5 7 14 2
o

1 6
Por ejemplo, el punto P2 = ( 5 , 5 ) se obtiene de resolver el siguiente sistema lineal:
ic


3x + 2y = 3
8x + 2y = 4
ér

que corresponden a los hiperplanos (rectas) que definen las restricciones (r1) y (r3).
Observación: el punto (0, 32 ) es solución de
um


3x + 2y = 3
x = 0
correspondiente a las restricciones (r1) y (r4) pero no es un vértice de la región factible Ω.
.N

Definición 4. Sea Ω una región poliedral cerrada de Rn determinada por un sistema de m

inecuaciones lineales. Se llaman vértices de Ω a los puntos pertenecientes a Ω que satisfa-
cen uno de los posibles sistemas de n ecuaciones lineales independientes (obtenido de las m
A

inecuaciones).

Ejercicio: determinar todos los vértices de la región poliedral cerrada Ω en R3 definida por


 x+y+z ≤ 3
y−z ≤ 2


 x − 2y ≤ 1
x ≥ 0


Geométricamente, los vértices de una región poliedral cerrada Ω de Rn son los puntos extre-
mos, esto es, puntos de Ω que no están contenidos en el interior de algún segmento contenido
en la región. (Formalizaremos la definición de punto extremo en la próxima clase).

5
 Clase 20 - Programación lineal (2)

El problema

Recordemos la formulación del problema de programación lineal

21
Minimizar f (x)
sujeto a Cx = d
Rx ≥ s

20
x≥0

o en su forma estándar
Minimizar cT x

F
sujeto a Ax = b
x≥0

A
con b ≥ 0.

M
En la clase anterior dijimos que la región factible es una región poliedral cerrada y que los
vértices de esa región serán importantes al momento de buscar las soluciones del problema
de PL.
A
-F

Método gráfico para el caso bidimensional

Cuando el problema tiene 2 variables la forma más simple de buscar la solución de un pro-
blema de PL es el método gráfico. Para fijar ideas, consideraremos el siguiente ejemplo en
o

dimensión 2:
book
ic

Minimizar z = −x1 − 2x2

2008/10
sujeto a −2x1 + x2 ≤ 2 page 98
−x1 + x2 ≤ 3
ér

x1 ≤ 3
x1 , x2 ≥ 0.
um

Se puede pensar en la función objetivo como z = f (x1 , x2 ) = cT (x1 , x2 ), con c = (−1, −2).
La región factible de este problema puede verse en la Figura 1.
.N

z = -15

6 - 2 x1 + x2 = 2

5
A

3
x1 =3
2
z = -4
- x1 + x2 = 3
1
z = -2

-3 -2 -1 0 1 2 3 4
-1 z =0

Graphical
Figura 1: Solución solution
gráfica of aproblema
de un linear program.
lineal.

and the parallel line 0 passes through the origin. The goal of the linear program is to 1
minimize the value of . As the figure illustrates, decreases as these lines move upward
and to the right. The objective cannot be decreased indefinitely, however. Eventually the
 Notar que la figura incluye algunas líneas punteadas correspondientes a las curvas de nivel
para diferentes valores de la función objetivo. Por ejemplo, el conjunto de nivel {(x1 , x2 )|z =
−x1 − 2x2 = −2} es la recta z = −2 que pasa por los puntos (2, 0) y (0, 1). El conjunto de ni-
vel {(x1 , x2 )|z = −x1 −2x2 = 0} es la recta paralela z = 0 que pasa por el origen. Recordemos
que el objetivo del problema es minimizar z. Como se ve en la Figura 1, el valor de z decrece
a medida que se avanza hacia la derecha, sin embargo no puede decrecer indefinidamente
porque la región factible es acotada. Al continuar trazando rectas, asociadas a conjuntos de

21
nivel, hacia la derecha llegará un momento donde esa recta intersecará por última vez al con-
junto factible. En este problema el mínimo ocurre cuando z = −15 en el punto (3, 6), que
corresponde a un vértice del conjunto factible. No es coincidencia que sea un vértice, como

20
veremos más adelante en algunos resultados.
Por otro lado se sabe, por un resultado de Análisis de varias variables, que la dirección
del gradiente de una función es la dirección de máximo crecimiento y, análogamente, la

F
dirección de menos gradiente es la dirección de máximo descenso. Esto será muy útil cuando
el objetivo sea maximizar o minimizar una función lineal. Ahora bien, como la función es

A
lineal, el gradiente no es más que el vector de costos c. En este caso, c = (−1, −2) y por lo
tanto, si queremos minimizar la función z = −x1 − 2x2 , basta con trazar rectas paralelas en

M
la dirección de −c, es decir, (1, 2) hasta intersecar por última vez a la región factible. Ver
Figura 2:
A
-F
o
ic
ér
um
.N

Figura 2: Solución gráfica de un problema lineal y vector c.

A

En resumen, el método gráfico para problemas de PL consiste en trazar rectas perpendicula-

res al vector gradiente de la función objetivo, es decir al vector de costos c, y trasladarse en
una dirección u otra dependiendo si se desea minimizar o maximizar.

Tipos de solución

Es fácil imaginar que pueden surgir diferentes tipos de problemas y soluciones para el caso
de un problema de PL en dos dimensiones, debido a la geometría del conjunto factible, a la
dirección del vector de costos y si se está minimizando o maximizando. En todas las figuras
siguientes consideraremos f (x1 , x2 ) = x1 + 2x2 .

2
 - Regiones no acotadas sin vértices:

21
20
Figura 3: Regiones no acotadas sin vértices.

F
En el caso de R1 se alcanza el valor mínimo en toda la recta (frontera de R1) y no hay valor

A
máximo. En el caso de R2 no hay mínimo ni máximo. Ver Figura 3.
- Regiones no acotadas con vértices:

M
A
-F
o
ic

Figura 4: Regiones no acotadas con vértices.

ér

En el caso de R3 no hay máximo ni mínimo. En el caso de R4 se alcanza el mínimo en los

vértices A y B, por lo tanto en todo el segmento que une a estos puntos, y no hay máximo.
um

Notar que en el caso de la región R4 hay infinitas soluciones al problema de minimización.

Ver Figura 4.
- Región acotada (con al menos tres vértices):
.N
A

Figura 5: Región acotada (con al menos tres vértices).

Asume el mínimo en el punto A y el máximo en el punto C de la región R5. Ver Figura 5.

3
 - Casos degenerados:

21
20
Figura 6: Casos degenerados

F
En el caso de la región R6 no hay mínimo ni máximo. En cambio, la región R7 tiene un

A
máximo en el punto A y no hay mínimo. Ver Figura 6.

M
A
-F
o

Figura 7: Casos degenerados

ic

La región R8 tiene un mínimo en el punto B y un máximo en el punto A. La región R9 tiene

ér

un mínimo y un máximo en el punto A. Ver Figura 7.

um

El caso n-dimensional

El método gráfico presentado en la sección anterior es útil sólo en casos bidimensionales

y en algunos casos tridimensionales. Para dimensiones mayores que 2 el método gráfico
.N

se torna impracticable y se requiere un método eficiente para resolver tales problemas. Pa-
ra esto vamos a introducir el Método Simplex, el cual fue propuesto en 1947 por George
Dantzig, para resolver problemas de PL que modelizaban situaciones de planificación eco-
A

nómica y militar (planning). Este método consiste en un algoritmo eficiente de búsqueda,

que comenzando en algún vértice de la región factible avanza hacia otro vértice vecino hasta
encontrar la solución óptima de una manera inteligente y eficiente. La Programación Lineal
se había desarrollado muy poco hasta 1947 debido a la dificultad computacional de resol-
ver problemas lineales por el hecho de tener una gran número de combinaciones posibles
de restricciones a ser consideradas. Con la aparición de este método, estos cálculos fueron
optimizados y desde entonces el método Simplex es uno de los algoritmos más estudiados y
utilizados a nivel mundial. Aunque se han propuesto otros métodos más recientes, el método
Simplex sigue siendo muy competitivo y eficiente y se puede aplicar en una gran cantidad
de problemas. En la actualidad, matemáticos aplicados, programadores y especialistas en
computación continuan investigando en mejores implementaciones y variantes del método.

4
 Clase 21 - Programación lineal (3)
Fundamentos matemáticos del método Simplex

Por razones de tiempo veremos algunos fundamentos matemáticos del método Simplex y
posteriormente presenteramentos la versión tabla del método, sin entrar en casos particulares

21
ni variantes especiales que lo hacen más robusto. Comenzaremos dando algunas definiciones.
Consideremos el sistema lineal de ecuaciones de la forma

20

 a11 x1 + a12 x2 + . . . + a1n xn = b1

 a21 x1 + a22 x2 + . . . + a2n xn = b2

.. .. .. .. .. .. .. .. .. (1)


 . . . . . . . . .

F
am1 x1 + am2 x2 + . . . + amn xn = bm


Matricialmente se puede expresar como

A
A x = b,

M
con A ∈ Rm×n , x ∈ Rn y b ∈ Rm . Supongamos que n > m y que la matriz A tiene rango m.
Definición 1. Solución básica: dado un conjunto de m ecuaciones lineales en n incógnitas
A
como en (1), se llamará base a cualquier submatriz B, m × m de A, no singular, formada
por m columnas independientes de A. Es posible reordenar columnas de A para definir B.
-F

Como rango(A) = m, es claro que existe tal matriz B. Sea N la submatriz m × (n − m) de A,

formada por las restantes columnas. Así podemos escribir A y x particionados:
 
xB
A = [B|N] y x = y por lo tanto, BxB + NxN = b.
o

xN
ic

Luego xB = B−1 b − B−1 NxN .

Si x = (xB , xN ) = (B−1 b − B−1 NxN , xN ) es tal que xN = 0, es decir cada una de las
ér

(n − m) componentes son iguales a 0, entonces xB = B−1 b y así x = (xB , xN ) = (xB , 0) es

llamada solución básica de (1), con respecto a la base B. Las componentes de xB son lla-
um

madas variables básicas de x y las componentes de xN son llamadas variables no-básicas

de x.
Definición 2. Solución básica degenerada: si una o más variables básicas en una solución
básica toma el valor 0, se dice que esta solución es una solución básica degenerada. Desde
.N

el punto de vista geométrico, esto ocurre cuando se tiene un vértice determinado por la
intersección de más de dos rectas o hiperplanos.
A

Ahora consideremos el siguiente sistema


Ax = b
(2)
x ≥ 0
con A ∈ Rm×n , de rango m y b ∈ Rm , el cual representa el conjunto factible de restricciones
de un problema de PL en su forma estándar.
Definición 3. Una solución x que satisface (2) se dice que es factible para este sistema. Si
además esa solución x es básica se llama solución básica factible. Si esta solución fuera
básica degenerada, se la llama solución básica factible degenerada.

1
 Finalmente, consideremos el problema de PL en la forma estándar:

 Minimizar cT x
sujeto a Ax = b (3)
x≥0


21
donde A ∈ Rm×n es una matriz de rango m.

Definición 4. Se llama solución básica factible óptima a la solución básica factible que da

20
el valor óptimo (en este caso el valor mínimo) para la función objetivo del problema (3).

El siguiente teorema, que sólo enunciaremos, establece la importancia fundamental de las

F
soluciones básicas factibles para resolver un problema de PL. Quien esté interesado en ver la
demostración puede consultarla en el libro “Linear and Nonlinear Programming” de Luen-

A
berger, Ye, 4th edition, Springer, 2015.

M
Teorema 1. Teorema Fundamental de Programación Lineal. Dado un problema de PL
en la forma estándar (3), donde A ∈ Rm×n , tiene rango m. Luego,
A
i) si existe una solución factible, entonces existe una solución básica factible.
-F

ii) si existe una solución factible óptima, entonces existe una solución básica factible
óptima.
o

La importancia de este teorema radica en que reduce el problema de resolver un problema

de PL a la búsqueda en el conjunto de soluciones básicas factibles. Como para un problema
ic

de PL de n variables y m restricciones, con n > m, se tiene (a lo sumo)

ér

 
n n!
=
m m!(n − m)!
um

soluciones básicas (correspondientes al número de posibilidades de elegir m columnas de las

n posibles para construir la matriz base B). De esta manera, este teorema da una natural, pe-
ro muy ineficiente, técnica de búsqueda finita, la cual suele llamarse método enumerativo.
Este método consiste en determinar todas las soluciones básicas factibles resolviendo todos
.N

los respectivos sistemas lineales y elegir el que produzca el menor valor objetivo, en el caso
de un problema de minimización. Por otro lado, el método simplex, que veremos próxima-
mente realiza una búsqueda inteligente sin necesidad de revisar todas las soluciones básicas
A

factibles.

Por último veremos la conexión de las soluciones básicas factibles con los puntos extremos.
Esta conexión establece la relación entre resultados algebraicos con resultados geométricos
y de convexidad.

Definición 5. Un punto x en la región factible Ω = {x|A x = b, x ≥ 0}, con A ∈ Rm×n de rango

m y b ∈ Rm , se dice que es un punto extremo de Ω si y sólo si x no puede ser expresado
como:
x = α y + (1 − α) z, y, z ∈ Ω, 0 < α < 1, x 6= y 6= z

2
 Esta definición dice básicamente que un punto extremo no puede estar en un segmento que
conecte otros dos puntos distintos del mismo conjunto. Por ejemplo, cada vértice de un
triángulo es un punto extremo del triángulo. Por lo tanto, todos los vértices de una región
poliedral cerrada como Ω son puntos extremos.
Ahora veremos algunos resultados que relacionan los vértices de la región factible con la
solución óptima de un problema de PL.

21
Teorema 2. Sea Ω la región factible dada en la definición anterior. Un vector x es un punto
extremo de Ω si y sólo es una solución básica factible de Ω.

20
Demostración.
(⇐) Supongamos que x es una solución básica factible de Ω, entonces x puede particionarse
como x = (xB , xN ), con BxB = b y xB ≥ 0, xN = 0.

F
Demostraremos que x debe ser un punto extremo, por contradicción. Si x no es un punto

A
extremo entonces existen dos puntos distintos y, z ∈ Ω tal que

M
x = αy + (1 − α)z, 0 < α < 1.

Como α y (1 − α) son positivos y xN = 0 entonces los vectores y, z también tienen se parti-

A
cionan de la misma manera y tienen n − m componentes iguales a cero:
-F

y = (yB , yN ), z = (zB , zN )

con ByB = b, yB ≥ 0, yN = 0, y BzB = b, zB ≥ 0, zN = 0.

Como la matriz B es inversible y Bx = By = Bz = b, entonces x = y = z lo cual contradice
o

que x no es punto extremo. Por lo tanto hemos probado que x es un punto extremo.
ic

(⇒) Supongamos que x es un punto extremo entonces veamos que x es una solución básica
factible. También se probará por contradicción.
ér

Como x es un punto extremo de Ω también debe ser un punto factible, es decir, Ax = b y

x ≥ 0. Reordenando las variables del vector x, si fuera necesario, es posible escribir x como
um

 
xB
x= ,
xN
.N

donde xN = 0 y xB > 0, es decir xN contiene las componentes de x que son iguales a cero.
Entonces particionamos la matriz A de igual forma como

A = [B|N]
A

donde B y N contienen las columnas correspondientes a xB y xN , respectivamente. Notar que

B podría no ser una matriz cuadrada. Si las columnas de B son linealmente independientes,
entonces x es una solución básica factible y termina aquí la demostración. Ahora supongamos
que las columnas de B son linealmente dependientes y construiremos dos puntos factibles
1 1
distintos y, z tal que x = y + z, y por lo tanto, x no puede ser un punto extremo.
2 2
Sea Bi la i-ésima columna de B. Si las columnas de B son linealmente dependientes, entonces
existen números reales r1 , . . . , rk , no todos nulos, tales que

r1 B1 + r2 B2 + · · · + rk Bk = 0, o equivalentemente B1 r1 + B2 r2 + · · · + Bk rk = 0,

3
 esto dice que si r = (r1 , r2 , . . . , rk ), entonces Br = 0. Notar que

B(xB ± αr) = BxB ± αBr = BxB ± 0 = BxB = b,

para todo α ∈ R. Como xB > 0, para valores suficientemente pequeños de ε > 0 se tiene que
xB + εr > 0 y xB − εr > 0.
Sean

21
   
xB + εr xB − εr
y= yz= .
xN xN
1 1

20
Luego, los vectores y, z son factibles, distintos entre sí y distintos de x. Como x = y + z,
2 2
lo cual contradice que x es un punto extremo.

Los corolarios que se enuncian a continuación son consecuencia de este teorema y del Teo-

F
rema Fundamental de PL.

A
Corolario 1. Si el conjunto Ω es no vacío, entonces tiene al menos un punto extremo.

M
Corolario 2. Si existe una solución óptima de un problema de programación lineal, entonces
existe una solución óptima finita que es un punto extremo del conjunto de restricciones.
Corolario 3. El conjunto de restricciones Ω tiene a lo sumo un número finito de puntos
A
extremos.
-F

Ejemplo: para fijar ideas vamos a determinar las soluciones básicas factibles del siguiente
problema de maximización.
Maximizar z = 4500x1 + 8000x2
o

sujeto a 5x1 + 20x2 ≤ 400 (r1)

10x1 + 15x2 ≤ 450
ic

(r2)
x1 ≥ 0 (r3)
x2 ≥ 0 (r4)
ér

Agregando variables de holgura, las restricciones de este problema en formato estándar re-
sultan:
um


 5x1 + 20x2 + x3 = 400
10x1 + 15x2 + x4 = 450
x1 , x2 , x3 , x4 ≥ 0


Así la matriz A y el vector b del conjunto de restricciones en el formato estándar están dados
.N

por    
5 20 1 0 400
A= b= .
10 15 0 1 450
A

 
4 4!
Como A tiene 2 filas y 4 columnas, se tienen a lo sumo = = 6 casos a
2 2!(4 − 2)!
considerar, como posibles vértices del conjunto factible. Ver Figura 1.
Caso 1: La matriz base B formada por las columnas 3 y 4 de la matriz A. Entonces
   
1 0 x3
B= xB = .
0 1 x4
Luego     
1 0 x3 400
B xB = b ⇒ =
0 1 x4 450

4
 21
20
F
Figura 1: Región factible del ejemplo.

A
Por lo tanto xB = (x3 , x4 ) = (400, 450), xN = (x1 , x2 ) = (0, 0), y la solución básica asociada

M
a esta base es (x1 , x2 , x3 , x4 ) = (0, 0, 400, 450), que corresponde al vértice A.
Caso 2: La matriz base B formada por las columnas 2 y 4 de la matriz A. Entonces
A
   
20 0 x2
B= xB = .
-F

15 1 x4

Luego     
20 0 x2 400
B xB = b ⇒ =
15 1 x4 450
o

Por lo tanto xB = (x2 , x4 ) = (20, 150), xN = (x1 , x3 ) = (0, 0), y la solución básica asociada a
ic

esta base es (x1 , x2 , x3 , x4 ) = (0, 20, 0, 150), que corresponde al vértice B.

ér

Caso 3: La matriz base B formada por las columnas 1 y 2 de la matriz A. Entonces

   
5 20 x1
um

B= xB = .
10 15 x2

Luego     
5 20 x1 400
B xB = b ⇒ =
.N

10 15 x2 450
Por lo tanto xB = (x1 , x2 ) = (24, 14), xN = (x3 , x4 ) = (0, 0), y la solución básica asociada a
esta base es (x1 , x2 , x3 , x4 ) = (24, 14, 0, 0), que corresponde al vértice C.
A

Caso 4: La matriz base B formada por las columnas 1 y 3 de la matriz A. Entonces

   
5 1 x1
B= xB = .
10 0 x3

Luego     
5 1 x1 400
B xB = b ⇒ =
10 0 x3 450
Por lo tanto xB = (x1 , x3 ) = (45, 175), xN = (x2 , x4 ) = (0, 0), y la solución básica asociada a
esta base es (x1 , x2 , x3 , x4 ) = (45, 0, 175, 0), que corresponde al vértice D.

5
 Caso 5: La matriz base B formada por las columnas 2 y 3 de la matriz A. Entonces
   
20 1 x2
B= xB = .
15 0 x3

Luego     
20 1 x2 400
B xB = b ⇒ =

21
15 0 x3 450
Por lo tanto xB = (x2 , x3 ) = (30, −200), xN = (x1 , x4 ) = (0, 0), y la solución básica asociada
a esta base es (x1 , x2 , x3 , x4 ) = (0, 30, −200, 0), que corresponde al punto E, que no es un

20
vértice.
Caso 6: La matriz base B formada por las columnas 1 y 4 de la matriz A. Entonces

F
   
5 0 x1
B= xB = .
10 1 x4

A
Luego     
B xB = b
M
⇒
5 0
10 1
x1
x4
=
400
450
Por lo tanto xB = (x1 , x4 ) = (80, −350), xN = (x2 , x3 ) = (0, 0), y la solución básica asociada
A
a esta base es (x1 , x2 , x3 , x4 ) = (80, 0, 0, −350), que corresponde al punto F, que no es un
vértice.
-F
o
ic
ér
um
.N
A

6
 Clase 22 - Programación lineal (4)
El método Simplex

El método Simplex es un método iterativo para resolver un problema de PL en la forma

estándar. Si el problema es no degenerado, el método va recorriendo la región factible, de

21
una manera eficiente, de una solución básica factible (vértice) a otra hasta llegar a la solución
básica factible óptima. En cada iteración, se verifica si la solución actual es óptima. Si no
lo es, el método elige una dirección donde se mejore el valor óptimo y se avanza en esa

20
dirección hasta encontrar otra solución factible óptima. Luego se repite este procedimiento.
Para entender mejor el método consideraremos un ejemplo:

minimizar z = −x1 − 2x2

F
sujeto a − 2x1 + x2 ≤ 2

A
−x1 + 2x2 ≤ 7
x1 ≤ 3

M
El formato estándar de este problema está dado por
x1 , x2 ≥ 0
A
minimizar z = −x1 − 2x2
-F

sujeto a − 2x1 + x2 + x3 = 2
−x1 + 2x2 + x4 = 7
x1 + x5 = 3
o

x1 , x2 , x3 , x4 , x5 ≥ 0.
ic

Para este problema se tienen:

ér

   
−2 1 1 0 0 2
A =  −1 2 0 1 0  y b =  7 .
1 0 0 0 1 3
um

Cuando se tienen restricciones del tipo ≤ y se agregan variables de holgura, se puede obtener
fácilmente una solución básica factible tomando xB = (x3 , x4 , x5 ) = (2, 7, 3) y xN = (x1 , x2 ) =
(0, 0), que corresponde al origen de coordenadas, es decir, el punto xa = (x1 , x2 , x3 , x4 , x5 ) =
.N

(0, 0, 2, 7, 3) en la Figura 1.
Para pasar a otra solución básica factible, alguna de las variables no básicas pasará a ser
básica, y por lo tanto alguna de las variables básicas dejará de serlo y se convertirá en no
A

básica. Para eso notemos que, a partir de la forma estándar, se obtiene que

x3 = 2 + 2 x1 − x2
x4 = 7 + x1 − 2 x2 (1)
x5 = 3 − x1
y
z = −x1 − 2 x2 ,
cuyo valor en xa es z = 0, pues x1 = x2 = 0. Como el objetivo es minimizar z, si cualquiera
de estas dos variables toma valores positivos z decrecerá. Es claro que conviene elegir la
variable x2 porque el coeficiente de x2 hará que z decrezca más que si se eligiera x1 . Por

1
 lo tanto, la variable x2 que era no básica (x2 = 0) pasará a ser básica, es decir tomará un
valor positivo. Debido a las restricciones (1), esta variable no puede crecer indefinidamente.
Manteniendo la variable no básica x1 , es decir x1 = 0, y teniendo en cuenta (1), obtenemos

x3 = 2 − x2 (2)
x4 = 7 − 2 x2 (3)

21
x5 = 3. (4)

Estas 3 ecuaciones permitirán determinar cual es la variable básica que pasará de ser básica book

20
a no básica: la variable x2 puede aumentar hasta que x3 o x4 se hagan cero sin tomar valores 2008/10/23
negativos. De (2), cuando x3 = 0, entonces x2 = 2 y además x4 = 3 > 0. De (3), cuando page 127
x4 = 0, entonces x2 = 7/2 y además x3 = −3/2 < 0. Por lo tanto, x2 puede crecer hasta el
valor x2 = 2, y x3 dejará de ser variable básica para ser no básica. Y así, la nueva solución

F
básica factible será el punto xb = (x1 , x2 , x3 , x4 , x5 ) = (0, 2, 0, 3, 3) en la Figura 1.

A
M x2
6

5
A
xd
xc
xf 4
-F

3
xb
2

1
o

xa xe
ic

-7 -6 -5 -4 -3 -2 -1 1 2 3 4 x1
ér

The simplex method.

Figura 1: El método Simplex.
We now test if this point is optimal. To do this, we determine if there exist any feasible
um

descent directions. The constraints can be written with the basic variables expressed in terms
Paraof the nonbasic la
terminar variables:
iteración, resta actualizar el problema en términos de la nueva solución
básica factible, donde las variables2 básicas
3 2 1 2están en xB = (x2 , x4 , x5 ) y las no básicas xN =
(x1 , x3 ). Para esto, se deben expresar 4 7 las 22
1 nuevas variables básicas en términos de las nuevas
3
.N

variables no básicas. Las variables 5 1

no básicas deben aparecer en la función objetivo y en la
matriz de restricciones se debe obtener una matriz
All other feasible points can be found by varying the values of the identidad en las variables
nonbasic variables 1 and
básicas.
and using the constraints to determine the values of the basic variables 3 , 4 , and 5 .
la2 primera
De Because ecuación
the nonbasic de (1),
variables tenemos
are currently zeroque = 2 + 2x
and xall2 variables 1 −be
must , y reemplazando
x3nonnegative, it en la función
is only yvalid
entolas
increase a nonbasic variable so that itde
becomes positive.
A

objetivo restricciones, el problema PL resulta

Our goal is to minimize the objective function
1 2 2 z = −4 − 5x1 + 2x3
minimizar
and its current value is 1 or 2 a
0. If eithersujeto is increased
x2 = from
2 + 2x 1−
zero, then
x3 will decrease.
This shows that nearby feasible points obtained by xincreasing either 1 or 2 give lower
4 = 3 − 3x1 + 2x3
values of the objective function, so that the current basis is not optimal.
The simplex method moves from one basis to an x5 adjacent
= 3 − basis,
x1 deleting and adding
just one variable from the basis. This corresponds to moving between adjacent basic feasible
solutions. In geometric terms, the simplexx1 , x2method , x5 ≥along
, x3 , x4moves 0. edges of the feasible region.
It is not difficult to calculate how the objective function changes along an edge of the feasible
region,
Como xNand
= this
(x1contributes
, x3 ) = (0,to 0),
the simplicity
entoncesof zthe=method.
−4 y x = (x , x4 , x5 ) = (2, 3, 3). Esto termina la
For this example, moving to an adjacent basic feasibleBsolution2 corresponds to in-
primera
creasingiteración
either 1 or del2 ,método.
but not both. The coefficient of 2 is greater in absolute value
than the coefficient of 1 , so decreases more rapidly when 2 is increased. In the hope of
making more rapid progress toward the solution, we choose to increase 2 rather than 1 . 2
 Fórmulas generales

Consideremos el problema de PL en el formato estándar

minimizar z = cT x
sujeto a Ax = b
x ≥ 0.

21
Si x = (xB , xN ) es una solución básica factible para alguna submatriz B de A = [B|N], entonces
z = cTB xB + cTN xN

20
BxB + NxN = b, ⇒ xB = B−1 b − B−1 NxN ,

entonces reemplazando en z, se obtiene z = cTB B−1 b + (cTN − cTB B−1 N)xN .

F
Ahora, el valor actual de las variables básicas y la función objetivo se obtienen tomando
xN = 0, y denotamos xB = b̂ = B−1 b, y ẑ = cTB B−1 b.

A
Si definimos cTB B−1 = yT , es decir y = (cTB B−1 )T = B−T cB , entonces

M z = yT b + (cTN − yT N)xN .

Se llama costo reducido de x j a ĉ j , la entrada j-ésima en el vector ĉTN ≡ (cTN − cTB B−1 N).
A
Luego, z = ẑ + ĉTN xN .
-F

Test de optimalidad: se debe analizar que ocurre con la función objetivo si se aumenta
cada variable no básica a partir del valor cero. Si ĉ j > 0 la función objetivo aumentará; si
ĉ j = 0 la función objetivo se mantiene constante; si ĉ j < 0 la función objetivo disminuirá si
x j comienza a crecer a partir de cero. Si la solución no es óptima, entonces se puede elegir
o

una variable no básica xt con ĉt < 0 para entrar a la base.

ic

¿Cuál variable dejará de ser básica?: para esto recordemos que xB = B−1 b − B−1 NxN , y
dado que las variables no básicas son iguales a cero, excepto xt , entonces
ér

xB = b̂ − Ât xt ,
donde Ât = B−1 At , y At es la columna t de A. Ahora, como hicimos en el ejemplo anterior,
um

miremos la i-ésima componente de xB :

(xB )i = b̂i − âi,t xt .
.N

Si âi,t > 0, entonces (xB )i decrecerá a medida que xt aumente, y (xB )i será igual a cero cuando
xt = b̂i /âi,t . Si âi,t < 0, entonces (xB )i aumentará y si âi,t = 0 entonces (xB )i se mantiene
constante. Se aplicará el siguiente test del cociente mínimo para decidir cual variable básica
pasará a ser no básica:
A

( )
b̂i
xs = mı́n | âi,t > 0 .
1≤i≤m âi,t

Para determinar los nuevos valores de las variables básicas y la función objetivo se debe
calcular:
xB ← xB − Ât xs , y ẑ ← ẑ + ĉt xs .

Si âi,t ≤ 0 para todo i entonces ninguna de las variables básicas decrecerá a medida que xt
aumenta, por lo tanto xt puede crecer tanto como se quiera. Esto significa, que la función
objetivo decrecerá mientras xt → ∞, o sea que no tendrá un mínimo finito, y por lo tanto el
problema es no acotado.

3
 El Algoritmo Simplex

Se comienza con una matriz base B correspondiente a una solución básica factible xB = b̂ =
B−1 b ≥ 0. Los siguientes tres pasos resumen el método Simplex.
Paso 1: test de optimalidad.
Calcular yT = cTB B−1 ; ĉTN = cTN − yT N.

21
Si ĉTN ≥ 0, entonces la solución es óptima. Sino, elegir una variable xt tal que ĉt < 0, para
entrar a la base (pasar a ser básica).

20
Paso 2: test del cociente mínimo.
Calcular Ât = B−1 At . Determinar un índice s tal que

F
( )
b̂s b̂i
= mı́n | âi,t > 0 .

A
âs,t 1≤i≤m âi,t

M
Este test determina cual es la variable que abandona la base.
Si âi,t ≤ 0 para todo i, el problema es no acotado.
A
Paso 3: actualización.
-F

Actualizar la matriz base B y las correspondientes variables básicas xB .

Volver al Paso 1.
o

Observación 1: esta es una versión simplificada del método simplex. Existen variantes al-
gorítmicas para evitar algunos problemas que pueden aparecer como ciclos o restricciones
ic

redundantes. Además, existen diferentes implementaciones del método simplex, tanto gra-
tuitas como comerciales.
ér

Observación 2: a veces no es posible disponer de una solución básica facible inicial fácil-
mente. Para esto existen algunos métodos que se usan previamente para obtener esta solución
um

inicial (método de las dos fases, método de la M grande).

Observación 3: bajo adecuadas hipótesis el algoritmo del método simplex converge a una
solución básica factible óptima o determina que el problema es no acotado.
.N

Observación 4: cuando las variables involucradas son enteras, el problema es mucho más
complejo, tanto matemática como computacionalmente, y esto se estudia en una subárea de
Optimización llamada Programación Lineal Entera.
A

Formulación tabla (tableau) del método Simplex

Para resolver problemas grandes de PL, las implementaciones computacionales eficientes

del método Simplex se basan en el algoritmo descripto arriba. En cambio, para problemas
pequeños es conveniente usar la formulación tabla por ser una manera compacta y sistemática
de organizar las cuentas del métoodo Simplex.

4
 Para la formulación estándar del ejemplo:

minimizar z = −x1 − 2x2

sujeto a − 2x1 + x2 + x3 = 2
−x1 + 2x2 + x4 = 7
x1 + x5 = 3
x1 , x2 , x3 , x4 , x5 ≥ 0.

21
La tabla inicial está dada por:

20
base x1 x2 x3 x4 x5 LD
−z −1 −2 0 0 0 0

F
x3 −2 1 1 0 0 2 (5)

A
x4 −1 2 0 1 0 7
x5 1 0 0 0 1 3

M
LD, significa lado derecho en la primera fila de la tabla. Notar que en la segunda fila escribi-
mos −z. Esto se debe a que igualamos a cero a la función objetivo:
A
-F

z = −x1 − 2x2 ⇒ −z − x1 − 2x2 + 0x3 + 0x4 + 0x5 = 0.

La primera columna indica las variables básicas y en las 3 últimas filas se indican los coefi-
o

cientes de las restricciones.

ic

La tabla del problema original y en la base actual de cada iteración están dadas por:
ér

base xB xN LD base xB xN LD
−z cB cN 0 −z 0 cTN − cTB B−1 N −cTB B−1 b
B−1 N B−1 b
um

xB B N b xB I

Como en la Tabla 5, los coeficientes de las variables no básicas en la primera fila (costos
.N

reducidos) son negativos, significa que esa solución no es óptima. Elegimos x2 por tener
el costo reducido de mayor magnitud (−2). Esa será la variable no básica que entrará a la
base. Para decidir cual variable básica dejará de serlo calculamos el mínimo de los siguientes
cocientes:
A

 
2 7
mı́n , = 2,
1 2

como este valor de mínimo corresponde a la variable básica x3 , esa variable abandonará la
base y será reemplazada por x2 .

5
El paso siguiente consiste en transformar la tabla de manera de obtener una matriz identi-
dad en las nuevas variables básicas x2 , x4 , x5 . Antes de reemplazar x3 se debe identificar la
columna de x2 en la tabla

⇓
base x1 x2 x3 x4 x5 LD
−z −1 −2 0 0 0 0
⇒ x3 −2 1 1 0 0 2
x4 −1 2 0 1 0 7
x5 1 0 0 0 1 3

21
 
0
 1 
y transformarla, usando operaciones elementales por filas, en   0  . En este caso, se debe


20
0
sumar 2 veces la fila de x3 a la fila de −z y restar 2 veces la fila de x3 a la fila de x4 . Así se
obtiene la siguiente tabla

F
A
base x1 x2 x3 x4 x5 LD
−z −5 0 2 0 0 4
x2
x4
x5
−2 1

1 0
1 0 0
3 0 −2 1 0
0 0 1
2
3
3
M
A
-F

Ahora se comienza la segunda iteración. En la primera fila, el costo reducido de la variable

x1 es −5 < 0, por lo tanto esta solución básica factible no es óptima y x1 puede pasar a
ser básica (entrar a la base). Aplicando el test del cociente mínimo, se deduce que x4 es la
variable que abandonará la base:
o
ic

⇓
base x1 x2 x3 x4 x5 LD
−z −5 0
ér

2 0 0 4
x2 −2 1 1 0 0 2
⇒ x4 3 0 −2 1 0 3
um

x5 1 0 0 0 1 3

Luego, repitiendo este procedimiento, se obtienen las siguientes tablas, de las iteraciones tres
y cuatro:
.N

⇓
base x1 x2 x3 x4 x5 LD
−z 0 0 −4/3 5/3 0 9
A

x2 0 1 −1/3 2/3 0 4
x1 1 0 −2/3 1/3 0 1
⇒ x5 0 0 2/3 −1/3 1 2

6
 Finalmente,
base x1 x2 x3 x4 x5 LD
−z 0 0 0 1 2 13
x2 0 1 0 1/2 1/2 5
x1 1 0 0 0 1 3
x3 0 0 1 −1/2 3/2 3

Vemos que en esta última tabla, los costos reducidos de las variables no básicas (x4 y x5 ) son
ambos positivos, y por lo tanto esta solución básica factible es óptima. La solución final está
claramente expresada en la última columna: z = −13, x2 = 5, x1 = 3 y x3 = 3, y las variables
no básicas x4 = x5 = 0, lo que corresponde al punto xd de la Figura 1.

21
20
F
A
M
A
-F
o
ic
ér
um
.N
A