Curso de Métodos Numéricos para

¿ingenieros?

Pedro Fortuny Ayuso

CURSO 2022/23, EPIG, GIJÓN. UNIVERSIDAD DE OVIEDO
Email address: [email protected]
⃝
CC⃝
BY: Copyright © 2023 Pedro Fortuny Ayuso

This work is licensed under the Creative Commons Attribution 3.0

License. To view a copy of this license, visit
http://creativecommons.org/licenses/by/3.0/es/
or send a letter to Creative Commons, 444 Castro Street, Suite 900,
Mountain View, California, 94041, USA.
 Índice general

Introducción 5
1. Algunos comentarios 5
Capítulo 1. Aritmética finita y análisis del error 7
1. Notación exponencial 7
2. Un ejemplo detallado de coma flotante 8
3. El error, definiciones básicas 15
4. [Contenido no esencial] Acotar el error 18
Capítulo 2. Ecuaciones no lineales: soluciones numéricas 21
1. Introducción 21
2. Algoritmos “geométricos” 22
3. El algoritmo de bisección 23
4. El algoritmo de Newton-Raphson 23
5. El algoritmo de la secante 25
6. [Contenido no esencial] Puntos fijos 27
7. Velocidad de convergencia de Newton-Raphson 30
8. Anexo: Código en Matlab/Octave 32
Capítulo 3. Solución aproximada de sistemas lineales 35
1. Algunos ejemplos importantes 35
2. El algoritmo de Gauss y la factorizacón LU 48
3. Algoritmos aproximados (de punto fijo) 61
4. Anexo: Código en Matlab/Octave 63
Capítulo 4. Interpolación 67
1. Interpolación lineal (a trozos) 67
2. ¿Tendría sentido interpolar con parábolas? 68
3. Splines cúbicos: curvatura continua 69
4. El polinomio interpolador de Lagrange: un solo polinomio
para todos los puntos 75
5. Interpolación aproximada 77
6. Código de algunos algoritmos 81
Capítulo 5. Derivación e Integración Numéricas 87
1. Derivación Numérica: un problema inestable 87
3
4 ÍNDICE GENERAL

2. Integración Numérica: fórmulas de cuadratura 89

Capítulo 6. Ecuaciones diferenciales 97
1. Introducción 97
2. Lo más básico 99
3. Discretización 100
4. [Contenido no esencial] Errores: runcamiento y redondeo 103
5. [Contenido no esencial] Solución de EDOs e integral 105
6. El método de Euler: integral usando el extremo izquierdo 105
7. [Contenido no esencial] Euler modificado: “punto medio”106
8. Método de Heun (regla del trapecio) 108
9. El modelo Predador-Presa (o de Lotka-Volterra) 108
 Introducción

Estas notas cubren esencialmente los contenidos de las clases de

Teoría de la asignatura de Métodos Numéricos de la EPI de Gijón,
al estilo del autor. Como toda colección informal de notas, pueden
contener errores: errores de concepto, erratas, e incluso errores en los
enunciados y demostraciones. Aun así, he ido tratando de revisarlas
cada curso y corregir las faltas que encontraba. De todos modos,
Estas notas no necesariamente suplen un buen libro (o dos)

1. Algunos comentarios
El objetivo que pretendo cuando explico esta asignatura es triple:
• Introducir los problemas básicos que afronta el Cálculo Nu-
mérico, en su versión más sencilla.
• Acostumbrar a los alumnos a enunciar algoritmos con pre-
cisión y a estudiar su complejidad.
• Transmitir a los alumnos la importancia de comprobar que la
salida de un programa tiene sentido. Esto es mucho más impor-
tante de lo que parece. Que “un ordenador diga algo” no
significa nada, de por sí.
Durante las notas haré referencia a varias herramientas informáticas.
De momento:
Matlab: Aunque procuraré que todo el código que aparezca
pueda utilizarse con Octave, haré referencia casi exclusiva-
mente a Matlab, que es el programa al que están acostum-
brados los alumnos en la Universidad de Oviedo y que usa-
rán presumiblemente en otras asignaturas.
Wolfram Alpha: Es necesario conocer esta herramienta, pues
permite gran cantidad de operaciones (matemáticas y no)
en una simple página web.
Excel: Bastantes de los algoritmos pueden implementarse en
hojas de cálculo. Los alumnos deberían intentar hacerlo: así
se darán cuenta de la potencia que tiene una herramienta tan
común (y que ya saben manejar).
5
6 INTRODUCCIÓN

Quizás lo que más me interesa es que los alumnos

• Aprendan a enunciar algoritmos con precisión, con la con-
ciencia de que un algoritmo debe ser claro, finito y debe termi-
nar.
• Sean capaces de seguir un algoritmo y de entender su fun-
cionamiento.
Por supuesto, habrá que entender los algoritmos que se expongan y
ser capaces de seguirlos paso a paso (esto es un requisito que estimo
imprescindible). Todos los que se explicarán serán de memorización
sencilla porque o bien serán muy breves (la mayoría) o bien tendrán
una expresión geométrica evidente (y por tanto, su expresión ver-
bal será fácilmente memorizable). La precisión formal es mucho más
importante en esta asignatura que las “ideas geométricas”, porque
precisamente esta asignatura estudia los métodos formales de resol-
ver ciertos problemas, no su expresión geométrica o intuitiva.
Las secciones ó párrafos indicados como [Contenido no esencial]
solo se incluyen como información extra: no es mi intención evaluar
al alumno sobre su contenido pero pienso que son lo suficientemente
interesantes como para que a alguien le sirvan para aprender un poco
más.
 CAPÍTULO 1

Aritmética finita y análisis del error

Se revisan muy rápidamente las notaciones exponenciales, la co-

ma flotante de doble precisión, la noción de error y algunas fuentes
comunes de errores en los cómputos.

1. Notación exponencial
Los numeros reales pueden tener un número finito o infinito de
cifras. Cuando se trabaja con ellos, siempre se han de aproximar a
una cantidad con un número finito (y habitualmente pequeño) de
dígitos. Además, conviene expresarlos de una manera uniforme para
que su magnitud y sus dígitos significativos sean fácilmente recono-
cibles a primera vista y para que el intercambio de datos entre má-
quinas sea determinista y eficiente. Esto ha dado lugar a la notación
conocida como científica o exponencial. En estas notas se explica de
una manera somera; se pretende más transmitir la idea que describir
con precisión las reglas de estas notaciones (pues tienen demasiados
casos particulares como para que compense una explicación). A lo
largo de este curso, utilizaremos las siguientes expresiones:
DEFINICIÓN 1. Una notación exponencial de un número real es una
expresión del tipo
± A.B × 10C
10C × ± A.B
± A.BeC
donde A, B y C son números naturales (que pueden ser nulos), ± in-
dica un signo (que puede omitirse si es +). Cualquiera de esas expre-
siones se refiere al número ( A + 0.B) · 10C (donde 0.B es el número
“cero coma B…”.
Por ejemplo:
• El número 3.123 es el propio 3.123.
• El número 0.01e − 7 es 0.000000001 (hay ocho ceros después
de la coma y antes del 1).
7
8 1. ARITMÉTICA FINITA Y ANÁLISIS DEL ERROR

• El número 103 × −2.3 es −2300.

• El número −23.783e − 1 es −2.3783.
Pero por lo general, la notación científica asume que el número a la
derecha del punto decimal es un solo dígito no nulo:
DEFINICIÓN 2. La notación exponencial estándar es la notación ex-
ponencial en la que A es un número entre 1 y 9.
Finalmente, las máquinas almacenan —generalmente— los nú-
meros de una manera muy concreta, en coma flotante:
DEFINICIÓN 3. Un formato en coma flotante es una especificación de
una notación exponencial en la que la longitud de A más la longitud
de B es una cantidad fija y en la que los exponentes (el número C)
varían en un rango determinado.
Por tanto, un formato en coma flotante solo puede expresar una
cantidad finita de números (aquellos que puedan escribirse según la
especificación del formato).

2. Un ejemplo detallado de coma flotante

Los ordenadores trabajan con números decimales utilzando (en
binario, pero para la explicación es mejor utilizar base 10) dos forma-
tos:
• Punto “fijo” (fixed point): los números se almacenan en una
lista de k e elementos para la parte entera y k d para la par-
te decimal, donde ambos números k e y k d están fijos de an-
temano: Las operaciones con este tipo de aritmética tienen

± e1 e2 e3 e4 e5 e6 e7 d1 d2 d3 d4

FIGURA 1. Punto fijo con 7 dígitos enteros y 4 decimales

(y uno de signo).

siempre la misma precisión “absoluta” (cada operación se

realiza exactamente con k d decimales). Los valores extremos
que pueden almacenarse son ±10ke y los valores más peque-
ños significativos son ±10−kd . Cuando un valor “externo” se
almacena en una variable de este tipo, los dígitos decimales
más allá del k d se truncan. Por tanto, se puede llegar a per-
der una precisión de ≃ 10−kd como mucho. El problema de
este tipo de aritmética es doble:
 2. UN EJEMPLO DETALLADO DE COMA FLOTANTE 9

(1) Los extremos máximos por lo general son insuficientes

(en binario, 32 dígitos solo cubren unos 4 × 109 valores,
cuatro mil millones, que en seguida se queda corte).
(2) La resolución decimal puede también quedarse corta
en problemas de muy alta precisión (que hoy en día
son habituales en física e ingenierías).
• “Punto flotante” (usualmente llamado coma flotante, en cas-
tellano), (floating point): los números se almacenan en una
pareja formada por k m elementos de la llamada “mantisa” y
k e elementos del “exponente” (aparte del signo). En princi-
pio la estructura es similar al punto fijo: Un número en esta

± m1 m2 m3 m4 m5 m6 m7 m8 ± e1 e2

FIGURA 2. Punto flotante con mantisa de 8 dígitos ente-

ros y exponente de 3 (dos más signo), y un signo.

expresión (la de la figura 2) es igual a

±0.m1 m2 · · · m8 × 10±e1 e2
donde el signo en cada parte será el que esté especificado. Se
exige que el primer dígito de la mantisa sea no nulo (m1 ̸=
0). En binario, el signo + es un 0 y el singo − es un 1. Por de-
finición, el valor máximo que se puede almacenar en esta co-
ma flotante es ±0.99999999 × 1099 (un valor enorme), mien-
tras que en decimal, los valores ínfimos son ±0.1 × 10−99 (un
cero seguido de 100 ceros y un 1 tras la coma). Como se ve,
este tipo de almacenamiento permite cálculos en un rango
mucho más amplio que el punto fijo. Esto hace que sea po-
sible utilizar este tipo de aritmética en problemas en que los
valores de las variables sean muy altos y en problemas en
que sean muy bajos, sin perder precisión relativa. Cada vez
que un número (real) se almacena en esa coma flotante, al
ser de la forma
0.m1 · · · m8 × 10e
el error que se comete es del orden de 0.00000001 × 10e , es
decir, 10e−8 : por tanto, si e es grande, el error absoluto es gran-
de, y si e es pequeño, el error absoluto es pequeño. Lo que
caracteriza la coma flotante es que al almacenar un número en
coma flotante, el error relativo que se comete es siempre del mismo
orden.
10 1. ARITMÉTICA FINITA Y ANÁLISIS DEL ERROR

Pero para ello necesitamos dos definiciones. Supongamos que x0

es un número que se quiere calcular (o almacenar) y que x0 es una
aproximación.
DEFINICIÓN 4. El error absoluto cometido al utilizar x0 en lugar de
x0 es: | x0 − x0 |. El error relativo es
| x0 − x 0 |
| x0 |
(que solo tiene sentido si x0 ̸= 0, claro).
El error absoluto mide la “distancia” que hay entre el valor y su
aproximación. El error relativo mide el tamaño del error “compara-
do con” el valor exacto. Por lo general, el error relativo es más infor-
mativo que el absoluto, pero es importante tener en cuenta ambos,
siempre.
Por tanto, en coma flotante, el error relativo cometido al truncar
un valor y almacenarlo es siempre menor que 10−km (en decimal, en
binario será con un 2): viene dado por la longitud de la mantisa. El
problema es que el error absoluto puede ser muy grande (del orden
de 1099−8 = 1097 , cuando se almacenan números de ese tamaño.
EJEMPLO 1 (Desbordamiento). Se tiene un controlador digital que
almacena números en punto fijo (decimal) con 5 dígitos enteros y
4 decimales. El contador comienza en 0 y se aumenta en 0.001 cada
vez (suma milésimas de segundo). ¿Cuánto tardará en ocurrir un
desbordamiento?
Como se almacena en punto fijo con 5 + 4 dígitos, el mayor va-
lor posible es 99999.9999. El desbordamiento se produce cuando el
contador llegue a 105 , por tanto, cuando pasen
105 /10−3 = 108
milésimas de segundo, es decir, 105 segundos, que son
105
≃ 1.157
24 · 60 · 60
días. Es decir, al cabo de 1 día y medio se producirá un error a ciencia
cierta.
EJEMPLO 2 (Pérdida de precisión). El mismo contador que antes
se almacena en un registro en coma flotante con 5 dígitos de mantisa
y 4 de exponente (contando el signo). ¿Cuándo se perderá precisión?
El número 0.001 en nuestra coma flotante es 1 × 10−3 . Para que las
milésimas sigan teniendo relevancia, el número mayor que se puede
 2. UN EJEMPLO DETALLADO DE COMA FLOTANTE 11

considerar que tenga cinco cifras es:

99.999 = 0.99999 × 102 .
Por tanto, cuando el contador llegue a 99.999, se producirá una pérdi-
da de precisión (en este caso total: ya no hay milésimas y el contador
se queda parado en 100). Es decir, al cabo de 100 segundos el conta-
dor falla totalmente.
El ejemplo 2 es extremo (el formato de coma flotante utilizado
es muy burdo porque la mantisa es muy pequeña) pero es lo que le
ocurrió al sistema Patriot en la Guerra del Golfo, esencialmente.
La precisión simple sigue existiendo: por un lado, las entidades fi-
nancieras tienen obligación de trabajar con un número específico de
decimales. Por otro lado, los microcontroladores están diseñados para
utilizar punto fijo porque la implementación de las operaciones bá-
sicas (y, sobre todo, de la suma y la resta) es mucho más sencilla (y,
por tanto, más económica).
El formato estándar de coma flotante es el documento conocido
como IEEE-754, que posee varias actualizaciones (el original es de
1987).
2.1. [Contenido no esencial] El formato binario de doble pre-
cisión IEEE-754. El formato de doble precisión es la especificación
del IEEE sobre la representación de números reales en secuencias de
16, 32, 64, 128 bits (y más), y su representación en formato decimal.
Se explican a continuación las propiedades principales de la doble
precisión en binario.
Para expresar números en doble precisión se utilizan 64 bits, es
decir, 64 dígitos binarios. El primero indica el signo (un 0 es signo
positivo, un 1, signo negativo). Los 11 siguientes se utilizan para re-
presentar el exponente como se indicará, y los 52 restantes se utilizan
para lo que se denomina la mantisa. De esta manera, un número en
doble precisión tiene tres partes: s, el signo (que será 0 ó 1), e, el expo-
nente (que variará entre 0 y 2047 (pues 211 = 2048), y m, un número
de 52 bits. Dados tres datos s, e, m, el número real N que representan
es:
• Si e ̸= 0 y e ̸= 2047 (si el exponente no es ningún valor
extremo), entonces
N = (−1)s × 2e−1023 × 1.m,
donde 1.m indica “uno-coma-m” en binario. Nótese, y esto
es lo importante, que el exponente no es el número representa-
do por los 11 bits de e, sino que “se desplaza hacia la derecha”.
12 1. ARITMÉTICA FINITA Y ANÁLISIS DEL ERROR

Un e = 01010101011, que en decimal es 683 representa real-

mente la potencia de 2 2683−1023 = 2−340 . Los e cuyo primer
bit es cero corresponden a potencias negativas de 2 y los que
tienen primer bit 1 a potencias positivas (210 = 1024).
• Si e = 0 entonces, si m ̸= 0 (si hay mantisa):
N = (−1)s × 2−1023 × 0.m,
donde 0.m indica “cero-coma-m” en binario.
• Los ceros con signo: si e = 0 y m = 0, el número es +0 o −0,
dependiendo de s. (es decir, el 0 tiene signo).
• El caso de e = 2047 (es decir, los 11 digitos del exponente
son 1) se reserva para codificar ±∞ y otros objetos que se
denominan NaN (Not-a-Number, que indica que una ope-
ración es ilegítima, como 1/0 o log(−2) o acos(3), en una
operación con números reales).
En realidad, el estándar es mucho más largo y completo, como es
natural, e incluye una gran colección de requisitos para los sistemas
electrónicos que realicen cálculos en coma flotante (por ejemplo, es-
pecifica cómo se han de truncar los resultados de las operaciones
fundamentales para asegurar que si se puede obtener un resultado
exacto, se obtenga).
Las ventajas de la coma flotante (y en concreto de la doble preci-
sión) son, aparte de su estandarización, que permite a la vez operar
con números muy pequeños (el número más pequeño que puede al-
macenar es 2−1023 ≃ 10−300 ) y números muy grandes (el mayor es
alrededor de 21023 ≃ 10300 ). La contrapartida es que, si se trabaja si-
multáneamente con ambos tipos de datos, los pequeños pierden pre-
cisión y desaparecen (se produce un error de cancelación o de trun-
cación). Pero si se tiene cuidado, es un formato enormemente útil y
versátil.

2.2. [Contenido no esencial] Conversión de base decimal a ba-

se dos y vuelta. Explico a continuación, para quien no sepa, cómo
se trasnforma un número en base decimal a base dos (binario) y al
revés.
En el pseudocódigo de estas notas, se utilizará la siguiente nota-
ción: la expresión x ← a significa que x es una variable, a representa
un valor (así que puede ser un número u otra variable) y a x se le
asigna el valor de a. La expresión u = c es la expresión condicional
¿es el valor designado por u igual al valor designado por c? Además,
m//n indica el cociente de dividir el número entero m ≥ 0 entre el
 2. UN EJEMPLO DETALLADO DE COMA FLOTANTE 13

número entero n > 0 y m%n es el resto de dicha divisón. Es decir,

m = (m//n) × n + (m%n).

Finalmente, si x es un número real x = A.B, la expresión { x } indica

la parte fraccionaria de x, es decir, 0.B.

Algoritmo 1 (Paso de decimal a binario, sin signo)

Input: A.B un número en decimal, con B finito, k un entero positi-
vo (que indica el número de digitos que se desea tras la coma en
binario)
Output: a.b (el número x en binario, truncado hasta 2−k )
if A = 0 then
a ← 0 y calular solo la PARTE DECIMAL
end if
⋆PARTE ENTERA
i ← −1, n ← A
while n > 0 do
i ← i+1
xi ← n%2
n ← n//2
end while
a ← xi xi−1 . . . x0 [la secuencia de restos en orden inverso]
⋆PARTE DECIMAL
if B = 0 then
b←0
return a.b
end if
i ← 0, n ← 0.B
while n > 0 and i < k do
i ← i+1
m ← 2n
if m ≥ 1 then
bi ← 1
else
bi ← 0
end if
n ← {m} [la parte decimal de m]
end while
b ← b1 b2 . . . bi
return a.b
14 1. ARITMÉTICA FINITA Y ANÁLISIS DEL ERROR

El Algoritmo 1 es una manera de pasar un número A.B en deci-

mal con un número finito de cifras a su forma binaria. El Algoritmo
2 se utiliza para realizar la operación inversa. Téngase en cuenta que,
puesto que hay números con una cantidad finita de digitos decimales
que no se pueden expresar con una cantidad finita de digitos en bina-
rio (el ejemplo más obvio es 0.1), se ha de especificar un número de
cifras decimales para la salida (así que no se obtiene necesariamente
el mismo número, sino una truncación).
Algoritmo 2 (Paso de binario a decimal, sin signo)
Input: A.B, un número positivo en binario, k un número entero no
negativo (el número de decimales que se quiere utilizar en bina-
rio).
Ouput: a.b, la truncación hasta precisión 2−k en decimal del núme-
ro A.B
⋆PARTE MAYOR QUE 0
Se escribe A = Ar Ar−1 . . . A0 (los dígitos binarios)
a ← 0, i ← 0
while i ≤ r do
a ← a + A i × 2i
i ← i+1
end while
⋆PARTE DECIMAL
if B = 0 then
return a.0
end if
b ← 0, i ← 0
while i ≤ k do
i ← i+1
b ← b + Bi × 2−i
end while
return a.b

Pasar de binario a decimal es “más sencillo” pero requiere ir su-

mando: no obtenemos un dígito por paso, sino que se han de sumar
potencias de 2. Se describe este proceso en el Algoritmo 2. Nótese que
el número de decimales de la salida no está especificado (se podría,
pero solo haría el algoritmo más complejo). Finalmente, en todos los
pasos en que se suma Ai × 2i o bien Bi × 2−i , tanto Ai como Bi son
o bien 0 o bien 1, así que ese producto solo significa “sumar o no
sumar” la potencia de 2 correspondiente (que es lo que expresa un
dígito binario, al fin y al cabo).
 3. EL ERROR, DEFINICIONES BÁSICAS 15

3. El error, definiciones básicas

Siempre que se opera con números con una cantidad finita de ci-
fras y siempre que se toman medidas en la realidad, hay que tener en
cuenta que contendrán, casi con certeza, un error. Esto no es grave.
Lo prioritario es tener una idea de su tamaño y saber que según vayan
haciéndose operaciones puede ir propagándose. Al final, lo que im-
porta es acotar el error absoluto, es decir, conocer un valor (la cota) que
sea mayor que el error cometido, para saber con certeza cuánto, como
mucho, dista el valor real del valor obtenido.
En lo que sigue, se parte de un valor exacto x (una constante, un
dato, la solución de un problema…) y de una aproximación, x̃.
DEFINICIÓN 5. Se llama error absoluto cometido al utilizar x̃ en lugar
de x al valor absoluto de la diferencia: | x − x̃ |.
Pero, salvo que x sea 0, uno está habitualmente más interesado en
el orden de magnitud del error, es decir, “cuánto se desvía x̃ de x en
proporción a x”:
DEFINICIÓN 6. Se llama error relativo cometido al utilizar x̃ en lugar
de x, siempre que x ̸= 0, al cociente
| x̃ − x |
|x|
(que es siempre un número positivo).
No vamos a utilizar una notación especial para ninguno de los
dos errores (hay autores que los llaman ∆ y δ, respectivamente, pero
cuanta menos notación innecesaria, mejor).
EJEMPLO 3. La constante π, que es la razón entre la longitud de la
circunferencia y su diámetro, es, aproximadamente 3.1415926534+
(con el + final se indica que es mayor que el número escrito hasta la
última cifra). Supongamos que se utiliza la aproximación π̃ = 3.14.
Se tiene:
• El error absoluto es |π − π̃ | = 0.0015926534+.
• El error relativo es |π −
π
π̃ |
≃ 10−4 × 5.069573.
Esto último significa que se produce un error de 5 diezmilésimas por
unidad (que es más o menos 1/2000) cada vez que se usa 3.14 en
lugar de π. Por tanto, si se suma 3.14 dos mil veces, el error cometido
al usar esa cantidad en lugar de 2000 × π es aproximadamente de 1.
Este es el significado interesante del error relativo: su inverso es el
número de veces que hay que sumar x̃ para que el error acumulado
sea 1 (que, posiblemente, será el orden de magnitud del problema).
16 1. ARITMÉTICA FINITA Y ANÁLISIS DEL ERROR

Antes de continuar analizando errores, conviene definir las dos

maneras más comunes de escribir números utilizando una cantidad
fija de dígitos: el truncamiento y el redondeo. No se van a dar las defi-
niciones precisas porque en este caso la precisión parece irrelevante.
Se parte de un número real (posiblemente con un número infinito de
cifras):
x = a1 a2 . . . ar .ar+1 . . . an . . .
Se definen:
DEFINICIÓN 7. El truncamiento de x a k cifras (significativas) es el
número a1 a2 . . . ak 0 . . . 0 (un número entero), si k ≤ r y si no, a1 . . . ar .ar+1 . . . ak
si k > r. Se trata de cortar las cifras de x y poner ceros si aun no se ha
llegado a la coma decimal.
DEFINICIÓN 8. El redondeo de x a k cifras (significativas) es el si-
guiente número:
• Si ak+1 < 5, entonces el redondeo es igual al truncamiento.
• Si 5 ≤ ak+1 ≤ 9, entonces el redondeo es igual al trunca-
miento más 10r−k+1 .
Este redondeo se denomina redondeo hacia más infinito, porque siem-
pre se obtiene un número mayor o igual que el truncamiento.
El problema con el redondeo es que pueden cambiar todas las
cifras. La gran ventaja es que el error que se comete al redondear es
menor que el que se comete al truncar (puede ser hasta de la mitad):
EJEMPLO 4. Si x = 178.299 y se van a usar 4 cifras, entonces el
truncamiento es x1 = 178.2, mientras que el redondeo es 178.3. El
error absoluto cometido en el primier caso es 0.099, mientas que en
el segundo es 0.001.
EJEMPLO 5. Si x = 999.995 y se van a usar 5 cifras, el truncamiento
es x1 = 999.99, mientras que el redondeo es 1000.0. Pese a que todas
las cifras son diferentes, el error cometido es el mismo (en este caso,
0.005) y esto es lo importante, no que los digitos “coincidan”.
¿Por qué se habla entonces de truncamiento? Porque cuando uno
trabaja en coma flotante, es inevitable que se produzca truncamiento
(pues el número de dígitos es finito) y se hace necesario tenerlo en
cuenta. Actualmente (2012) la mayoría de programas que trabajan
en doble precisión, realmente lo hacen con muchos más dígitos in-
ternamente y posiblemente al reducir a doble precisión redondeen.
Aun así, los truncamientos se producen en algún momento (cuando
se sobrepasa la capacidad de la máquina).
 3. EL ERROR, DEFINICIONES BÁSICAS 17

3.1. [Contenido no esencial] Fuentes del error. El error se origi-

na de diversas maneras. Por una parte, cualquier medición está sujeta
a él (por eso los aparatos de medida se venden con un margen esti-
mado de precisión); esto es intrínseco a la naturaleza y lo único que
puede hacerse es tenerlo en cuenta y saber su magnitud (conocer una
buena cota). Por otro lado, las operaciones realizadas en aritmética
finita dan lugar tanto a la propagación de los errores como a la apari-
ción de nuevos, precisamente por la cantidad limitada de dígitos que
se pueden usar.
Se pueden enumerar, por ejemplo, las siguientes fuentes de error:
• El error en la medición, del que ya se ha hablado. Es inevitable.
• El error de truncamiento: ocurre cuando un número (dato o
resultado de una operación) tiene más dígitos que los utili-
zados en las operaciones y se “olvida” una parte.
• El error de redondeo: ocurre cuando, por la razón que sea, se
redondea un número a una precisión determinada.
• El error de cancelación: se produce cuando una operación da
lugar a errores relativos mucho más grandes que los absolu-
tos. Habitualmente tiene lugar al sustraer de magnitud muy
similar. El ejemplo canónico de esto aparece en la inestabili-
dad de la ecuación cuadrática.
• El error de acumulación: se produce al acumular (sumar, esen-
cialmente) pequeños errores del mismo signo muchas veces.
Es lo que ocurrió en el suceso de los Patriot en febrero de
1991, en la operación Tormenta del Desierto1.
En la aritmética de precisión finita, todos estos errores ocurren. Con-
viene quizás conocer las siguientes reglas (que son los peores casos
posibles):
• Al sumar números del mismo signo, el error absoluto puede
ser la suma de los errores absolutos y el error relativo, lo
mismo.
• Al sumar números de signo contrario, el error absoluto se com-
porta como en el caso anterior, pero el error relativo puede au-
mentar de manera incontrolada: 1000.2 − 1000.1 solo tiene una
cifra significativa (así que el error relativo puede ser de has-
ta un 10%, mientras que en los operandos, el error relativo
era de 1 × 10−4 .

1Pueden consultarse el siguiente documento de la Univ. de Texas:

http://www.cs.utexas.edu/~downing/papers/PatriotA1992.pdf
y el informe oficial: http://www.fas.org/spp/starwars/gao/im92026.htm
18 1. ARITMÉTICA FINITA Y ANÁLISIS DEL ERROR

• Al multiplicar, el error absoluto tiene la magnitud del ma-

yor factor por el error en el otro factor. Si los factores son de
la misma magnitud, puede ser el doble del mayor error ab-
soluto por el mayor factor. El error relativo es de la misma
magnitud que el mayor error relativo (y si son de la misma
magnitud, su suma).
• Al dividir por un número mayor o igual que 1, el error absoluto
es aproximadamente el error absoluto del numerador parti-
do por el denominador y el error relativo es el error relativo
del numerador (esencialmente como en la multiplicación).
Al dividir por números cercanos a cero, lo que ocurre es que
se perderá precisión absoluta y si luego se opera con núme-
ros de magnitud similar, el error de cancelación puede ser
importante. Compárense las siguientes cuentas:
33
26493 − = −0.256 (el resultado buscado) .
0.0012456
33
26493 − = −13.024
0.001245
Si bien el error relativo de truncación es solo de 4 × 10−4
(media milésima), el error relativo del resultado es de 49.8
(es decir, el resultado obtenido es 50 veces más grande que
el que debía ser). Este (esencialmente) es el problema fun-
damental de los métodos de resolución de sistemas que uti-
lizan divisiones (como el de Gauss y por supuseto, el de Cra-
mer). Cuando se explique el método de Gauss, se verá que
es conveniente buscar la manera de hacer las divisiones con
los divisores más grandes posibles (estrategias de pivotaje).

4. [Contenido no esencial] Acotar el error

Como se dijo antes, lo importante no es saber con precisión el
error cometido en una medición o al resolver un problema, pues con
casi certeza, esto será imposible, sino tener una idea y, sobre todo, te-
ner una buena cota: saber que el error absoluto cometido es menor que
una cantidad y que esta cantidad sea lo más ajustada posible (sería
inútil, por ejemplo, decir que 2.71 es una aproximación de e con un
error menor que 400).
Así pues, lo único que se podrá hacer realmente será estimar un
número mayor que el error cometido y razonablemente pequeño. Eso
es acotar.
 4. [Contenido no esencial] ACOTAR EL ERROR 19

4.1. Algunas cotas. Lo primero que ha de hacerse es acotar el

error si se conoce una cantidad con una variación aproximada. Es
decir, si se sabe que x = a ± ϵ, donde ϵ > 0, el error que se comete al
utilizar a en lugar de x es desconocido (esto es importante) pero es a
lo sumo, ϵ. Por lo tanto, el error relativo que se comete es, a lo sumo el
error absoluto dividido entre el menor de los valores posibles en valor ab-
soluto: téngase en cuenta que esto es delicado, pues si a − ϵ < 0 pero
a + ϵ > 0 entonces no se puede acotar el error relativo porque x podría
ser 0 (y, como ya se explicó, el error relativo solo tiene sentido para
cantidades no nulas).
EJEMPLO 6. Si se sabe que π = 3.14 ± 0.01, se sabe que el error
absoluto máximo cometido es 0.01, mientras que el error relativo es
como mucho
0.01
≃ .003194
|3.13|
(más o menos 1/300).
Téngase en cuenta que, si lo que se busca es una cota superior y
se ha de realizar una división, ha de tomarse el divisor lo más pequeño
posible (cuanto menor es el divisor, mayor es el resultado). Con esto
se ha de ser muy cuidadoso.
Las reglas de la sección 3.1 son esenciales si se quiere acotar el
error de una serie de operaciones aritméticas. Como se dice allí, ha de
tenerse un cuidado muy especial cuando se realizan divisiones con
números menores que uno, pues posiblemente se llegue a resultados
inútiles (como que “el resultado es 7 pero el error absoluto es 23”).
 CAPÍTULO 2

Ecuaciones no lineales: soluciones numéricas

Se estudia en este capítulo el problema de resolver una ecuación

no lineal de la forma f ( x ) = 0, dada la función f y ciertas condiciones
inciales.
Cada uno de los algoritmos que se estudia tiene sus ventajas e
inconvenientes; no hay que desechar ninguno a priori simplemente
porque “sea muy lento” —esta es la tentación fácil al estudiar la con-
vergencia de la bisección. Como se verá, el “mejor” algoritmo de los
que se estudiarán —el de Newton-Raphson— tiene dos pegas impor-
tantes: puede converger (o incluso no hacerlo) lejos de la condición
inicial (y por tanto volverse inútil si lo que se busca es una raíz “cer-
ca” de un punto) y además, requiere computar el valor de la deriva-
da de f en cada iteración, lo cual puede tener un coste computacional
excesivo.
Antes de proseguir, nótese que todos los algoritmos se describen
con condiciones de parada en función del valor de f , no con condicio-
nes de parada tipo Cauchy. Esto se hace para simplificar la exposición.
Si se quiere garantizar un número de decimales de precisión, es ne-
cesario o bien evaluar cambios de signo cerca de cada punto o bien
hacer un estudio muy detallado de la derivada, que está más allá del
alcance de este curso a mi entender.

1. Introducción
Calcular raíces de funciones —y especialmente de polinomios—
es uno de los problemas clásicos de las Matemáticas. Hasta hace unos
doscientos años se pensaba que cualquier polinomio podía “resolver-
se algebraicamente”, es decir: dada una ecuación polinómica an x n +
· · · + a1 x + a0 = 0 habría una fórmula “con radicales” que daría to-
dos los ceros de ese polinomio, como la fórmula de la ecuación de se-
gundo grado pero para grado cualquiera. La Teoría de Galois mostró
que esta idea no era más que un bonito sueño: el polinomio general
de quinto grado no admite una solución en función de radicales.
21
22 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

De todos modos, buscar una fórmula cerrada para resolver una

ecuaión no es más que una manera de postponer el problema de la
aproximación, pues al fin y al cabo (y esto es importante):
Las únicas operaciones que se pueden hacer con precisión total
siempre son la suma, la resta y la multiplicación.
Ni siquiera es posible dividir y obtener resultados exactos1.
En fin, es lógico buscar maneras lo más precisas y rápidas posibles
de resolver ecuaciones no lineales y que no requieran operaciones
mucho más complejas que la propia ecuación (claro).
Se distinguiran en este capítulo dos tipos de algoritmos: los “geo-
métricos”, por lo que se entienden aquellos que se basan en una idea
geométrica simple y el “de punto fijo”, que requiere desarrollar la
idea de contractividad.
Antes de proseguir, hay dos requerimientos esenciales en cual-
quier algoritmo de búsqueda de raíces:
• Especificar una precisión, es decir, un ϵ > 0 tal que si | f (c)| <
ϵ entonces se asume que c es una raíz. Esto se debe a que, al
utilizar un número finito de cifras, es posible que nunca ocu-
rra que f (c) = 0, siendo c el valor que el algoritmo obtiene
en cada paso.
• Incluir una condición de parada. Por la misma razón o por fa-
llos del algoritmo o porque la función no tiene raíces, podría
ocurrir que nunca se diera que | f (c)| < ϵ. Es imprescindible
especificar un momento en que el algoritmo se detiene, de
manera determinista. En estas notas la condición de parada se-
rá siempre un número de repeticiones del algoritmo, N > 0:
si se sobrepasa, se devuelve un “error”.

2. Algoritmos “geométricos”
Si suponemos que la función f cuya raíz se quiere calcular es con-
tinua en un intervalo compacto [ a, b] (lo cual es una suposición razo-
nable), el Teorema de Bolzano puede ser de utilidad, si es que sus
hipótesis se cumplen. Recordemos:
TEOREMA (Bolzano). Sea f : [ a, b] → R una función continua en
[ a, b] y tal que f ( a) f (b) < 0 (es decir, cambia de signo entre a y b). En-
tonces existe c ∈ ( a, b) tal que f (c) = 0.

1Podría argumentarse que si se dividen números racionales se puede escri-

bir el resultado como decimales periódicos, pero esto sigue siendo “postponer el
problema”.
 4. EL ALGORITMO DE NEWTON-RAPHSON 23

El enunciado afirma que si f cambia de signo en un intervalo ce-

rrado, entonces al menos hay una raíz. Se podría pensar en mues-
trear el intervalo (por ejemplo en anchuras (b − a)/10i ) e ir buscan-
do subintervalos más pequeños en los que sigue cambiando el signo,
hasta llegar a un punto en que la precisión hace que hayamos locali-
zado una raíz “o casi”. En realidad, es más simple ir subdividiendo
en mitades. Esto es el algoritmo de bisección.

3. El algoritmo de bisección
Se parte de una función f , continua en un intervalo [ a, b], y de los
valores a y b. Como siempre, se ha de especificar una precisión ϵ > 0
(si | f (c)| < ϵ entonces c es una “raíz aproximada”) y un número
máximo de ejecuciones del bucle principal N > 0. Con todos estos
datos, se procede utilizando el Teorema de Bolzano:
Si f ( a) f (b) < 0, puesto que f es continua en [ a, b], tendrá alguna
a+b
raíz. Tómese c como el punto medio de [ a, b], es decir, c = . Hay
2
tres posibilidades:
• O bien | f (c)| < ϵ, en cuyo caso se termina y se devuelve el
valor c como raíz aproximada.
• O bien f ( a) f (c) < 0, en cuyo caso se sustituye b por c y se
repite todo el proceso.
• O bien f (b) f (c) < 0, en cuyo caso se sustituye a por c y se
repite todo el proceso.
La iteración se realiza como mucho N veces (para lo cual hay que
llevar la cuenta, obviamente). Si al cabo de N veces no se ha obtenido
una raíz, se termina con un mensaje de error.
El enunciado formal es el Algoritmo 3. Nótese que se utiliza un
signo nuevo: a ↔ b, que indica que los valores de a y b se intercam-
bian.

4. El algoritmo de Newton-Raphson
Una idea geométrica clásica (y de teoría clásica de la aproxima-
ción) es, en lugar de calcular una solución de una ecuación f ( x ) = 0
directamente, utilizar la mejor aproximación lineal a f , que es la recta
tangente en un punto. Es decir, en lugar de calcular una raíz de f ,
utilizar un valor ( x, f ( x )) para trazar la recta tangente a la gráfica de
f en ese punto y resolver la ecuación dada por el corte de esta tan-
gente con el eje OX, en lugar de f ( x ) = 0. Obviamente, el resultado
no será una raíz de f , pero en condiciones generales, uno espera que se
aproxime algo (la solución de una ecuación aproximada debería ser
24 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

Algoritmo 3 (Algoritmo de Bisección)

Input: una función f ( x ), un par de números reales, a, b con
f ( a) f (b) < 0, una tolerancia ϵ > 0 y un límite de iteraciones N > 0

Output o bien un mensaje de error o bien un número real r entre a

y b tal que | f (c)| < ϵ (una raíz aproximada)
⋆INICIO
r ← a+2 b
i←0
while | f (r )| ≥ ϵ and i < N do
if f (r ) f (b) < 0 then
a ← r [intervalo [r, b]]
else
b ← r [intervalo [ a, r ]]
end if
i ← i+1
r ← a+2 b
end while
if i ≥ N then
return ERROR
end if
return r

una solución aproximada). Si se repite el proceso de aproximación

mediante la tangente, uno espera ir acercándose a una raíz de f . Esta
es la idea del algoritmo de Newton-Raphson.
Recuérdese que la ecuación de la recta que pasa por el punto ( x0 , y0 )
y tiene pendiente b es:

Y = b ( X − x0 ) + y0

así que la ecuación de la recta tangente a la gráfica de f en el punto

( x0 , f ( x0 )) es (suponiendo que f es derivable en x0 ):

Y = f ′ ( x0 )( X − x0 ) + f ( x0 ).

El punto de corte de esta recta con el eje OX es, obviamente,

( )
f ( x0 )
( x1 , y1 ) = x0 − ′ ,0
f ( x0 )

suponiendo que existe (es decir, que f ′ ( x0 ) ̸= 0).

 5. EL ALGORITMO DE LA SECANTE 25

La iteración, si se supone que se ha calculado xn es, por tanto:

f ( xn )
x n +1 = x n − .
f ′ ( xn )
Así que se tiene el Algoritmo 4. En el enunciado solo se especifica un
posible lugar de error para no cargarlo, pero hay que tener en cuenta
que cada vez que se evalúa f o cada vez que se realiza una operación
(cualquiera) podría ocurrir un error de coma flotante. Aunque en los
enunciados que aparezcan de ahora en adelante no se mencionen to-
dos, el alumno debe ser consciente de que cualquier implementación
ha de emitir una excepción (raise an exception) en caso de que haya
un fallo de coma flotante.

Algoritmo 4 (Algoritmo de Newton-Raphson)

Input: una función f ( x ) derivable, una semilla x0 ∈ R, una toleran-
cia ϵ > 0 y un límite de iteraciones N > 0
Ouput: o bien un mensaje de error o bien un número real c tal que
| f (c)| < ϵ (es decir, una raíz aproximada)
⋆INICIO
i←0
while | f ( xi )| ≥ ϵ and i ≤ N do
f (x )
xi+1 ← xi − ′ i [posibles NaN e ∞]
f ( xi )
i ← i+1
end while
if i > N then
return ERROR
end if
return c ← xi

5. El algoritmo de la secante
f (x )
El algoritmo de Newton-Raphson contiene la evaluación f ′ ( xn ) ,
n
para la que hay que calcular el valor de dos expresiones: f ( xn ) y
f ′ ( xn ), lo cual puede ser exageradamente costoso. Además, hay mu-
chos casos en los que ni siquiera se tiene información real sobre f ′ ( x ),
así que puede ser hasta utópico pensar que el algoritmo es utilizable.
Una solución a esta pega es aproximar el cálculo de la derivada
utilizando la idea geométrica de que la tangente es el límite de las secan-
tes; en lugar de calcular la recta tangente se utiliza una aproximación
26 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

con dos puntos que se suponen “cercanos”. Como todo el mundo

sabe, la derivada de f ( x ) en el punto c es, si existe, el límite
f (c + h) − f (c)
lim .
h →0 h
En la situación del algoritmo de Newton-Raphson, si en lugar de uti-
lizar un solo dato xn , se recurre a los dos anteriores, xn y xn−1 , se
puede pensar que se está aproximando la derivada f ′ ( xn ) por medio
del cociente
f ( x n ) − f ( x n −1 )
f ′ ( xn ) ≃ ,
x n − x n −1
y por tanto, la fórmula para calcular el término xn+1 quedaría
x n − x n −1
x n +1 = x n − f ( x n ) ,
f ( x n ) − f ( x n −1 )
y con esto se puede ya enunciar el algoritmo correspondiente. Ténga-
se en cuenta que, para comenzarlo, hacen falta dos semillas, no basta
con una. El factor de f ( xn ) en la iteración es justamente el inverso de
la aproximación de la derivada en xn utilizando el punto xn−1 como
xn + h.

Algoritmo 5 (Algoritmo de la secante)

Input: Una función f ( x ), una tolerancia ϵ > 0, un límite de itera-
ciones N > 0 y dos semillas x−1 , x0 ∈ R
Output: O bien un número c ∈ R tal que | f (c)| < ϵ o bien un
mensaje de error
⋆INICIO
i←0
while | f ( xi )| ≥ ϵ and i ≤ N do
x i − x i −1
x i +1 ← x i − f ( x i ) [posibles NaN e ∞]
f ( x i ) − f ( x i −1 )
i ← i+1
end while
if i > N then
return ERROR
end if
return c ← xi

Es conveniente, al implementar este algoritmo, conservar en me-

moria no solo xn y xn−1 , sino los valores calculados (que se han uti-
lizado ya) f ( xn ) y f ( xn−1 ), para no computarlos más de una vez.
 6. [Contenido no esencial] PUNTOS FIJOS 27

6. [Contenido no esencial] Puntos fijos

Los algoritmos de punto fijo (que, como veremos, engloban a los
anteriores de una manera indirecta) se basan en la noción de contrac-
tividad, que no refleja más que la idea de que una función puede hacer
que las imágenes de dos puntos cualesquiera estén más cerca que los
puntos originales (es decir, la función contrae el conjunto inicial). Esta
noción, ligada directamente a la derivada (si es que la función es de-
rivable), lleva de manera directa a la de punto fijo de una iteración y,
mediante un pequeño artificio, a poder resolver ecuaciones generales
utilizando iteraciones de una misma función.
6.1. Contractividad: las ecuaciones g( x ) = x. Sea g una función
real de una variable real derivable en un punto c. Esto significa que,
dado cualquier infinitésimo o, existe otro o1 tal que
g(c + o) = g(c) + g′ (c)o + oo1 ,
es decir, que cerca de c, la función g se aproxima muy bien mediante una
función lineal.
Supongamos ahora que o es la anchura de un intervalo “peque-
ño” centrado en c. Si olvidamos el error supralineal (es decir, el tér-
mino oo1 ), pues es “infinitamente más pequeño” que todo lo demás,
se puede pensar que el intervalo (c − o, c + o) se transforma en el in-
tervalo ( g(c) − g′ (c)o, g(c) + g′ (c)o): esto es, un intervalo de radio o
se transforma aproximadamente, por la aplicación g en uno de an-
chura g′ (c)o (se dilata o se contrae un factor g′ (c)). Esta es la noción
equivalente de derivada que se utiliza en la fórmula de integración
por cambio de variable (el Teorema del Jacobiano): la derivada (y
en varias variables el determinante jacobiano) mide la dilatación o con-
tracción que sufre la recta real (un abierto de Rn ) justo en un entorno
infinitesimal de un punto.
Por tanto, si | g′ (c)| < 1, la recta se contrae cerca de c.
Supóngase que eso ocurre en todo un intervalo [ a, b]. Para no car-
gar la notación, supongamos que g′ es continua en [ a, b] y que se tiene
| g′ ( x )| < 1 para todo x ∈ [ a, b] —es decir, g siempre contrae la recta.
Para empezar, hay un cierto λ < 1 tal que | g( x )| < λ, por el Teorema
de Weierstrass (hemos supuesto g′ continua). Por el Teorema del Va-
lor Medio (en la forma de desigualdad, pero importa poco), resulta
que, para cualesquiera x1 , x2 ∈ [ a, b], se tiene que
| g( x1 ) − g( x2 )| ≤ λ| x1 − x2 |,
en castellano: la distancia entre las imágenes de dos puntos es menor que
λ por la distancia enre los puntos, pero como λ es menor que 1, resulta
28 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

que la distancia entre las imágenes es siempre menor que entre los puntos
iniciales. Es decir: la aplicación g está contrayendo el intervalo [ a, b] por
todas partes.
Hay muchos énfasis en el párrafo anterior, pero son necesarios
porque son la clave para el resultado central: la existencia de un punto
fijo.
Se tiene, por tanto, que la anchura del conjunto g([ a, b]) es menor
o igual que λ(b − a). Si ocurriera además que g([ a, b]) ⊂ [ a, b], es de-
cir, si g transformara el segmento [ a, b] en una parte suya, entonces
se podría calcular también g( g([ a, b])), que sería una parte propia
de g([ a, b]) y que tendría anchura menor o igual que λλ(b − a) =
λ2 (b − a). Ahora podría iterarse la composición con g y, al cabo de
n iteraciones se tendría que la imagen estaría contenida en un inter-
valo de anchura λn (b − a). Como λ < 1, resulta que estas anchuras
van haciéndose cada vez más pequeñas. Una sencilla aplicación del
Principio de los Intervalos Encajados probaría que en el límite, los
conjuntos g([ a, b]), g( g([ a, b])), …, gn ([ a, b]), …, terminan cortándo-
se en un solo punto α. Además, este punto, por construcción, debe
cumplir que g(α) = α, es decir, es un punto fijo. Hemos mostrado (no
demostrado) el siguiente
TEOREMA 1. Sea g : [ a, b] → [ a, b] una aplicación de un intervalo ce-
rrado en sí mismo, continua y derivable en [ a, b]. Si existe λ < 1 positivo tal
que para todo x ∈ [ a, b] se tiene que | g′ ( x )| ≤ λ, entonces existe un único
punto α ∈ [ a, b] tal que g(α) = α.
Así que, si se tiene una función de un intervalo en sí mismo cu-
ya deriva se acota por una constante menor que uno, la ecuación
g( x ) = x tiene una única solución en dicho intervalo. Además, la
explicación previa al enunicado muestra que para calcular α basta
con tomar cualquier x del intervalo e ir calculando g( x ), g( g( x )), …:
el límite es α, independientemente de x.
Por tanto, resolver ecuaciones de la forma g( x ) = x para funciones
contractivas es tremendamente simple: basta con iterar g.

6.2. Aplicación a ecuaciones cualesquiera f ( x ) = 0. Pero por lo

general, nadie se encuentra una ecuación del tipo g( x ) = x; lo que se
busca es resolver ecuaciones como f ( x ) = 0.
Esto no presenta ningún problema, pues:
f (x) = 0 ⇔ x − f (x) = x
de manera que buscar un cero de f ( x ) equivale a buscar un punto fijo de
g ( x ) = x − f ( x ).
 6. [Contenido no esencial] PUNTOS FIJOS 29

En realidad, si ϕ( x ) es una función que no se anula nunca, buscar

un cero de f ( x ) equivale a buscar un punto fijo de g( x ) = x − ϕ( x ) f ( x ).
Esto permite, por ejemplo, escalar f para que su derivada sea cercana
a 1 en la raíz y así que g′ sea bastante pequeño, para acelerar el proce-
so de convergencia. O se puede simplemente tomar g( x ) = x − c f ( x )
para cierto c que haga que la derivada de g sea relativamente pequeña
en valor absoluto.

6.3. El algoritmo. La implementación de un algoritmo de bús-

queda de punto fijo es muy sencilla. Como todos, requiere una to-
lerancia ϵ y un número máximo de iteraciones N. La pega es que el
algoritmo es inútil si la función g no envía un intervalo [ a, b] en sí
mismo. Esto es algo que hay que comprobar a priori:
NOTA 1. Sea g : [ a, b] → R una aplicación. Si se quiere buscar un
punto fijo de g en [ a, b] utilizando la propiedad de contractividad, es
necesario:
• Comprobar que g([ a, b]) ⊂ [ a, b].
• Si g es derivable2 en todo [ a, b], comprobar que existe λ ∈ R
tal que 0 < λ < 1 y para el que | g′ ( x )| ≤ λ para todo x ∈
[ a, b].
Supuesta esta comprobación, el algoritmo para buscar el punto
fijo de g : [ a, b] → [ a, b] es el siguiente:

6.4. Utilización para cómputo de raíces. El algoritmo de pun-

to fijo se puede utilizar para computar raíces, como se explicó en la
Sección 6.2; para ello, si la ecuación tiene la forma f ( x ) = 0, puede
tomarse cualquier función g( x ) de la forma
g( x ) = x − k f ( x )
donde k ∈ R es un número conveniente. Se hace esta operación para
conseguir que g tenga derivada cercana a 0 cerca de la raíz, y para
que, si χ es la raíz (desconocida), conseguir así que g defina una fun-
ción contractiva de un intervalo de la forma [χ − ρ, χ + ρ] (que será
el [ a, b] que se utilice).
Lo difícil, como antes, es probar que g es contractiva y envía un
intervalo en un intervalo.

2Si g no es derivable, se requiere una condición mucho más complicada de

verificar: que | g( x ) − g( x ′ )| ≤ λ| g( x ) − g( x ′ )| para cierto 0 < λ < 1 y para todo
par x, x ′ ∈ [ a, b]. Esto es la condición de Lipschitz con constante menor que 1.
30 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

Algoritmo 6 (Búsqueda de punto fijo)

Input: una función g (que ya se supone contractiva, etc…), una
semilla x0 ∈ [ a, b], una tolerancia ϵ > 0 y un número máximo de
iteraciones N > 0
Output: o bien un número c ∈ [ a, b] tal que |c − g(c)| < ϵ o bien
un mensaje de error
⋆INICIO
i ← 0, c ← x0
while |c − g(c)| ≥ ϵ and i ≤ N do
c ← g(c)
i ← i+1
end while
if i > N then
return ERROR
end if
return c

6.5. Velocidad de convergencia. Si se tiene la derivada acotada

en valor absoluto por una constante λ < 1, es relativamente sencillo
acotar la velocidad de convergencia del algoritmo del punto fijo, pues
como se dijo antes, la imagen del intervalo [ a, b] es un intervalo de
anchura λ(b − a). Tras iterar i veces, la imagen de la composición
gi ([ a, b]) = g( g(. . . ( g([ a, b])))) (una composición repetida i veces)
está incluida en un intervalo de longitud λi (b − a), de manera que se
puede asegurar que la distancia entre la composición gi ( x ) y el punto
fijo c es como mucho λi (b − a). Si λ es suficientemente pequeño, la
convergencia puede ser muy rápida.
EJEMPLO 7. Si [ a, b] = [0, 1] y |(| g′ ( x )) < .1, se puede asegurar que
tras cada iteración, hay un decimal exacto en el punto calculado, sea
cual sea la semilla.

7. Velocidad de convergencia de Newton-Raphson

Obsérvese que la expresión para calcular el punto siguiente en el
algoritmo de Newton-Raphson es:
f ( xn )
x n +1 = x n −
f ′ ( xn )
que corresponde a buscar un punto fijo de la función
1
g( x ) = x − ′ f ( x ).
f (x)
 7. VELOCIDAD DE CONVERGENCIA DE NEWTON-RAPHSON 31

que, como se explica arriba, es una manera de convertir una ecuación

f ( x ) = 0 en un problema de punto fijo. La derivada de g es, en este
caso:
f ′ ( x )2 − f ( x ) f ′′ ( x )
g′ ( x ) = 1 −
f ′ ( x )2
que, en el punto χ en que f (χ) = 0 vale g′ (χ) = 0. Es decir, en el
punto fijo, la derivada de g es 0. Esto hace que la convergencia (una
vez que uno está “cerca” del punto fijo) sea muy rápida (se dice de
segundo orden).
De hecho, se puede probar el siguiente resultado:
TEOREMA 2. Supongamos que f es una función derivable dos veces en
un intervalo [r − ϵ, r + ϵ], que r es una raíz de f y que se sabe que
• La derivada segunda está acotada superiormente: | f ′′ ( x )| < K
para x ∈ [r − ϵ, r + ϵ],
• La derivada primera está acotada inferiormente: | f ′ ( x )| > L para
x ∈ [r − ϵ, r + ϵ]
Entonces, si xk ∈ [r − ϵ, r + ϵ], el siguiente elemento de la iteración de
Newton-Raphson está también en el intervalo y
K
| r − x k +1 | < | ||r − xk |2 .
2L
Como corolario, se tiene que:
COROLARIO 1. Si las condiciones del Teorema 2 se cumplen y ϵ < 0.1 y
K < 2L, entonces a partir de k, cada iteración de Newton-Raphson obtiene
una aproximación con el doble de cifras exactas que la aproximación anterior.
PRUEBA. Esto es porque, si suponemos que k = 0, entonces x0 tiene
al menos una cifra exacta. Por el Teorema, x1 difiere de r en menos
que 0.12 = .01. Y a partir de ahora el número de ceros en la expresión
se va duplicando. □
Para utilizar el Teorema 2 o su corolario, es preciso:
• Saber que se está cerca de una raíz. La manera común de
verificarlo es computando f cerca de la aproximación (hacia
la derecha o la izquierda): si el signo cambia, ha de haber
una raíz (pues f es continua, al ser derivable).
• Acotar la anchura del intervalo (saber que se está a distancia
menor de 1/10, por ejemplo.
• Acotar f ′′ ( x ) superiormente y f ′ ( x ) inferiormente (esto pue-
de ser sencillo o no).
32 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

8. Anexo: Código en Matlab/Octave

Se incluyen a continuación algunos listados con implementacio-
nes “correctas” de los algoritmos descritos en este capítulo, utiliza-
bles tanto en Matlab como en Octave. Téngase en cuenta, sin embar-
go, que, puesto que ambos programas tienen cuidado de que las ope-
raciones en coma flotante no terminen en un error irrecuperable, las
implementaciones que se presentan no incluyen ninguna revisión de
posibles errores. Por ejemplo, en cualquiera de estos algoritmos, si se
utiliza la función log( x ) y se evalúa en un número negativo, el al-
goritmo posiblemente continuará y quizás incluso alcance una raíz
(real o compleja). No se han incluido estos tests para no complicar la
exposición.

8.1. Implementación del algoritmo de bisección. El siguiente có-

digo implementa el algoritmo de bisección en Octave/Matlab. Los
parámetros de entrada son:
f: una función anónima,
a: el extremo izquierdo del intervalo,
b: el extremo derecho del intervalo,
epsilon: una tolerancia (por defecto, eps),
n: un número máximo de iteraciones (por defecto, 50).
La salida puede ser nula (si no hay cambio de signo, con un mensaje)
o una raíz aproximada en la tolerancia o bien un mensaje (warning)
junto con el último valor calculado por el algoritmo (que no será una
raíz aproximada en la tolerancia). El formato de la salida, para faci-
litar el estudio es un par [z, N] donde z es la raíz aproximada (o el
valor más cercano calculado) y N es el número de iteraciones hasta
calcularla.

8.2. Implementación del algoritmo de Newton-Raphson. El al-

goritmo de Newton-Raphson es más sencillo de escribir (siempre sin
tener en cuenta las posibles excepciones de coma flotante), pero re-
quiere un dato más complejo en el input: la derivada de f , que debe
ser otra función anónima.
Como se ha de utilizar la derivada de f y no se quiere suponer
que el usuario tiene un programa con cálculo simbólico, se requiere
que uno de los parámetros sea explícitamente la función f ′ ( x ). En un
entorno con cálculo simbólico esto no sería necesario.
Así, la entrada es:
f: una función anónima,
fp: otra función anónima, la derivada de f,
 8. ANEXO: CÓDIGO EN MATLAB/OCTAVE 33

1 % Algoritmo de biseccion con error admitido y limite de parada

% Tengase en cuenta que en TODAS las evaluaciones f(...) podria
% ocurrir un error , esto no se tiene en cuenta en esta implementacion
% (en cualquier caso , se terimara)
function [z, N] = Bisec(f, x, y, epsilon = eps , n = 50)
N = 0;
if(f(x)*f(y) >0)
warning('no hay cambio de signo ')
return
end
11 % guardar valores en memoria
fx = f(x);
fy = f(y);
if(fx == 0)
z = x;
return
end
if(fy == 0)
z = y;
return
21 end
z = (x+y)/2;
fz = f(z);
while(abs(fz) >= epsilon & N < n)
N = N + 1;
% multiplicar SIGNOS , no valores
if(sign(fz)*sign(fx) < 0)
y = z;
fy = fz;
else
31 x = z;
fx = fz;
end
% podria haber un error
z = (x+y)/2;
fz = f(z);
end
if(N >= n)
warning('No se ha alcanzado el error admitido antes del limite.')
end
41 end

FIGURA 1. Código del algoritmo de Bisección

x0: la semilla,
epsilon: una tolerancia (por defecto eps),
N: el número máximo de iteraciones (por defecto 50).
El formato de la salida, para facilitar el estudio es un par [xn, N]
donde xn es la raíz aproximada (o el valor más cercano calculado) y
N es el número de iteraciones hasta calcularla.
 34 2. ECUACIONES NO LINEALES: SOLUCIONES NUMÉRICAS

% Implementacion del metodo de Newton -Raphson

function [z n] = NewtonF(f, fp , x0 , epsilon = eps , N = 50)
n = 0;
xn = x0;
% Se supone que f y fp son funciones
fn = f(xn);
while(abs(fn) >= epsilon & n <= N)
n = n + 1;
9 fn = f(xn); % evaluar una sola vez
% siguiente punto
xn = xn - fn/fp(xn); % podria haber una excepcion
end
z = xn;
if(n == N)
warning('No converje en MAX iteraciones ');
end
end

FIGURA 2. Implementación de Newton-Raphson

 CAPÍTULO 3

Solución aproximada de sistemas lineales

Se estudia en este capítulo una colección de algoritmos clásicos

para resolver de manera aproximada sistemas de ecuaciones linea-
les. Se comienza con el algoritmo de Gauss (y su interpretación como
factorización LU) y se discute brevemente su complejidad. Se apro-
vecha esta sección para introducir la noción de número de condición de
una matriz y su relación con el error relativo (de manera sencilla, sin
contemplar posibles errores en la matriz). A continuación se intro-
ducen los algoritmos de tipo punto fijo (se estudian específicamente
dos de ellos, el de Jacobi y el de Gauss-Seidel), y se compara su com-
plejidad con la del algoritmo de Gauss.
Aunque los algoritmos son interesantes de por sí (y su enuncia-
do preciso también), es igual de importante conocer las condiciones
de utilidad de cada uno, de manera al menos teórica (por ejemplo,
la necesidad de que la matriz de un algoritmo de punto fijo sea con-
vergente, noción análoga en este contexto a la contractividad). No se
discutirán en profundidad aspectos relacionados con el espectro de
una aplicación lineal, aunque es la manera más simple de compren-
der la necesidad de la condición de convergencia.
En todo este tema, se supone que ha de resolverse un sistema de
ecuaciones lineales de la forma

(1) Ax = b

donde A es una matriz cuadrada de orden n y b es un vector de Rn .

Se supone, siempre que A no es singular (es decir, que el sistema es
compatible determinado).

1. Algunos ejemplos importantes

Antes de proceder al estudio de cómo se resuelven los sistemas
de ecuaciones lineales numéricamente, conviene tener algún ejem-
plo importante y conocido. E, incluso antes de plantearse “resolver”
sistemas lineales, conviene conocer aplicaciones lineales útiles en la
vida actual.
35
36 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

A día de hoy, una de las aplicaciones más importantes del Ál-

gebra Lineal es el tratamiento de imágenes (en general, de señales).
Lo primero que vamos a estudiar es cómo puede transformarse una
imagen en un vector. El código necesario para las instrucciones de
Matlab que se indican está en mi página web.

1.1. Imágenes en “escala de grises”. Una imagen digital de n ×

m pixeles en color se compone (en el sistema RGB) de tres imágenes
monocromas, una para el color rojo, otra para el verde (green) y otra
para el azul (blue). En el modelo aditivo de color (que es el que utili-
zan las pantallas), con esos tres colores se pueden obtener (en teoría)
todos los del espectro.
Para simplificar el estudio, nosotros trabajaremos siempre con imá-
genes monocromáticas, que entenderemos en “escala de grises”.
La gran mayoría de las transformaciones de imágenes se tradu-
cen en realizar combinaciones lineales de los tonos de un punto con
los de otros. Incluso girar una imagen puede entenderse así. Pero pa-
ra trabajar cómodamente con estas transformaciones, es conveniente
convertir la imagen original en un vector: transformar la matriz de
tamaño n × m en un vector de nm coordenadas. Esto se realiza “con-
catenando” las filas de la imagen.
Llamemos I a la matriz que representa la imagen, un rectángulo
de n × m entradas en cada una de las cuales se almacena un valor
(que, de momento, supondremos que es un número entre 0 y 1). Así,
pues
I = ( Ijk ), j = 0, . . . , n − 1, k = 0, . . . , m − 1.
El primer subíndice indica la fila (hay n de ellas) y, el segundo, la
columna (hay m). Concatenemos las filas en un vector

v I = (vi ), i = 1, . . . , nm

de nm componentes. Concatenar las filas de I significa que: las pri-

meros n componentes de v I son los píxeles de la fila superior de I (se
suelen enumerar así). Los siguientes n corresponden a la fila siguien-
te (la que tiene j = 1), etc. hasta terminar. Formalmente:

vr = Ijk ⇔ r − 1 = jm + k.

Lo mejor para aclararse es hacer un ejemplo.

EJEMPLO 8. Consideremos una imagen de 4 filas por 6 columnas

como sigue (se han coloreado los primeros elementos de cada fila
 1. ALGUNOS EJEMPLOS IMPORTANTES 37

para que después se vea dónde caen en el vector):

1 2 0 4 6 8
2 3 1 6 5 4
6 4 3 2 2 1
9 2 3 1 4 0
Compruébese cómo, la entrada (1, 2) tiene el valor 1 (recuérdese que
se indexan los elementos desde 0). El vector correspondiente, de 24
componentes, es:
1 2 0 4 6 8 2 3 1 6 5 4 6 4 3 2 2 1 9 2 3 1 4 0
Como puede comprobarse, la componente 9 = 1 · 6 + 3 (es decir, v9 ,
contando como que v1 es la primera) es 1, como corresponde al píxel
de índices j = 1, k = 2, pues 9 − 1 = 1 · 6 + 2.
Una imagen de mayor tamaño solo da lugar a un vector más gran-
de pero no hay complicación niguna.
La transformación es bidireccional: un vector de nm componentes
se puede convertir en una matriz n × m haciendo. Para ello, si los
elementos de la matriz son Ijk , entonces Ijk = vr para r − 1 = jm + k,
donde k varía entre 0 y m − 1.
En todo este capítulo utilizaremos la fotografía que se muestra en
la Figura 1, que es una imagen de 512 filas y 331 columnas en escala de
grises de un ave (obtenida de wallhere.com/es/wallpaper/853988).
Para cargarla en Matlab, basta ejecutar

FIGURA 1. Fotografía utilizada en este capítulo, de ta-

maño 512 × 331.

>> load('ave.dat','-mat');
que crea una variable bird en la que están almacenados los valores
de la escala de grises de la imagen. Para visualizarla,
38 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

>> image(bird);
>> colormap(gray);
1.1.1. Convoluciones. Una de las técnicas más utilizadas en proce-
samiento de la imagen es la convolución. Este término tiene un signifi-
cado matemático preciso que no vamos a explicar, pues simplemente
lo aplicaremos a este problema gráfico.
En procesamiento de imágenes, se entiende por convlución la apli-
cación de un “núcleo”, sobre todos los puntos de la imagen. Un nú-
cleo es una matriz cuadrada “pequeña” N de número impar de filas
(habitualmente de tamaño 3 × 3 ó 5 × 5, aunque puede ser algo ma-
yor) que se “desplaza” sobre toda la imagen. Esta matriz N actúa
como un filtro y transforma la entrada Iij en la suma de los valores
de las entradas que rodean al Iij ponderadas por los valores de Nkl ,
como en la figura 2 (en realidad es un vídeo: en algún dispositivo
puede verse como tal).

FIGURA 2. Un núcleo sobre una imagen.

En dicha figura, la imagen I corresponde a la matriz “azul”, gran-

de, de tamaño 5 × 5, el núcleo es la matriz que aparece como “subín-
dices” y se desplaza, de tamaño 3 × 3; y la imagen obtenida J es más
pequeña que la inicial, de tamaño 3 × 3, en verde. Es decir:
 
3 3 2 1 0    
0 0 1 3 1 0 1 2 12 12 17
 
I=   
3 1 2 2 3 , K = 2 2 0 , J = 10 17 19
 
2 0 0 2 2 0 1 2 9 6 14
2 0 0 0 1
Por supuesto, un núcleo como el K de la imagen no tiene mucho senti-
do porque puede convertir números entre 0 y 1 en números mayores
que 1. Por lo general, los núcleo que se utilizan en procesamiento de
imagen tienen como peculiaridad que la suma de todas sus entradas
es como mucho 1: es decir, transforman el valor de un punto de la
imagen en una “media” ponderada de dicho punto y los que lo ro-
dean (los 8 de alrededor, en este caso).
 1. ALGUNOS EJEMPLOS IMPORTANTES 39

function [J] = nucleo1(I, K)

[f, c] = size(I);
3 J = zeros(size(I));
Ks = size(K,1)*size(K,2); % filas , columnas
Kv = reshape(K, 1, Ks);
f2 = (size(K,1) -1)/2; % filas
c2 = (size(K,2) -1)/2; % columnas
for k = f2+1:f-f2
for l = c2+1:c-c2
J(k,l) = Kv * reshape(I(k-1:k+1,l-1:l+1),Ks ,1);
end
end
13 J;
end

FIGURA 3. Una función para aplicar un núcleo a una

imagen (convolución).

1.1.2. Desenfoque. El ejemplo más sencillo de convolución es des-

enfocar una imagen. Puede entenderse que “desenfocar” consiste en
hacer que la intensidad de color de cada punto se redistribuya entre
él y los de alrededor. Por ejemplo, un núcleo como
 
0.025 0.05 0.025
K =  0.05 0.7 0.05 
0.025 0.05 0.025

haría que cada punto de la imagen perdiera el 30% de su intensidad y

la redistribuyera entre los de alrededor, dando el 5% a los de su norte,
sur, este y oeste, y el 2.5% a los de las esquinas. Comprobemos cómo
afecta este núcleo a la imagen de la Figura 1. Como vamos a iterar el
desenfoque, llamamos J a la imagen de partida y de final. Utilizamos
el programa nucleo1.m de la Figura 1.1.2 para aplicar el núcleo K a
la imagen J.
>> J = double(bird); % para poder multiplicar
>> K = [0.025 0.05 0.025; 0.05 0.7 0.05; 0.025 0.05 0.025];
>> J = nucleo1(J, K);
>> image(J);
Si solo se realiza una vez el desenfoque, casi no se percibe. Si se repite
el bucle unas cuantas veces, va observándose cómo se pierde claridad
en los contornos (que es, justamente, en lo que consiste desenfocar).
Este desenfoque es conocido como “gaussiano”, pues se parece a una
distribución “campana de Gauss” alrededor de cada punto.
 40 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

function [J] = nucleo(I, K)

[f, c] = size(I); % filas y columnas de una vez
[fn , cn] = size(K);
fn2 = (fn -1) /2;
cn2 = (cn -1) /2;
6 I1 = [zeros(fn2 ,c+cn -1); zeros(f,cn2) I zeros(f,cn2); zeros(fn2 ,c+cn
-1)];
J = zeros(size(I));
Ks = size(K,1)*size(K,2);
Kv = reshape(K, 1, Ks);
for k = fn2 +1:f+fn2;
for l = cn2 +1:c+cn2;
J(k-fn2 ,l-fn2) = Kv * reshape(I1(k-fn2:k+fn2 ,l-cn2:l+cn2),Ks
,1);
end
end
J;
16 end

FIGURA 4. Una función para aplicar un núcleo a una

imagen (convolución) ampliando la original con 0.

Otro desenfoque distinto consiste en “redistribuir” la intensidad

en un punto entre él y sus adyacentes homogéneamente:
 
1 1 1
1
K = 1 1 1
9 1 1 1

1.1.3. Los bordes de la imagen. Quien haya experimentado con las

operaciones anteriores y repetido varias veces la convolución, obser-
vará que la imagen modificada tiene bordes negros cada vez más an-
chos. Esto se debe a que, en la función nucleo1, los bucles no cubren
toda la fila y columna de J, sino solo la parte en la que el núcleo cabe
completamente. En la Figura 2 ya podía observarse este fenómeno y
cómo la imagen obtenida, tenía tamaño menor que la original.
Este problema se “resuelve” expandiendo la matriz original para
que el núcleo quepa en todos los puntos (por tanto, hay que añadir
columnas a la izquierda y a la derecha y filas arriba y abajo, tantas co-
mo la mitad del tamaño del núcleo menos uno). Nosotros haremos la
expansión más simple: por medio de ceros. Se implementa en el pro-
grama nucleo.m de la Figura 4. Como se ve en el listado, hace falta
ajustar bien el tamaño de la matriz sobre la que se aplica la convolu-
ción y, al final, devolver solo la parte del resultado que interesa. Este
ajuste es realmente importante en todos los métodos de convolución.
Nuestra solución es la más elemental (no necesariamente la peor).
 1. ALGUNOS EJEMPLOS IMPORTANTES 41

1.1.4. Expresión matricial general de las convoluciones. Lo interesan-

te de estas transformaciones es que el valor de cada elemento del
destino (cada punto de la imagen final) es una combinación lineal de
valores de elementos del origen. Por este motivo, la transformación
de convolución es una aplicación lineal: envía un vector (en este caso,
de tamaño nm) en otro del mismo tamaño (o menor si no se hace el
ajuste).

EJEMPLO 9. En nuestro caso, la imagen es demasiado grande pa-

ra escribir explícitamente la aplicación con una matriz. Tomemos un
caso más pequeño: una imagen de 5 × 4 píxeles. El espacio vectorial
al que pertenece tiene dimensión 20, así que la matriz de transforma-
ción será una matriz 20 × 20. Lo que vamos a ver es que pese a ser muy
grande, es muy simple.
Si realizamos la convolución con un núcleo 3 × 3:
 
k11 k12 k13
K = k21 k22 k23 
k31 k32 k33

para los elementos del interior de la imagen, Iij , se tiene, como se ve

en el listado de la Figura 4:
 
Ii−1,j−1
 Ii−1,j 
 
 .. 
Jij = (k11 k12 · · · k32 k33 )  . 
 
 Ii+1,j 
Ii+1,j+1

donde el vector fila son los elementos de K en filas y el vector columna

son los 9 (en este caso) puntos que circundan el Iij . Como se ve, el
elemento Jij solo está afectado por sí mismo y los otros 8 elementos
de I que “lo rodean”: es decir, en la matriz M de la transformación,
si v I es el vector correspondiente a la imagen I, cada fila solo tendrá
9 entradas: distribuidas según corresponde. Recuérdese que, si v I =
(vi ), entonces
vr = Iij ⇔ r − 1 = im + j
de donde la submatriz 3 × 3:
 
Ii−1,j−1 Ii−1,j Ii−1,j+1
 Ii,j−1 Iij Ii,j+1 
Ii+1,j−1 Ii+1,j Ii+1,j+1
42 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

se corresponde con las componentes siguientes del vector v I :

   
Ii−1,j−1 Ii−1,j Ii−1,j+1 v ( i −1) m + j v ( i −1) m + j +1 v ( i −1) m + j +2
 Ii,j−1 Iij Ii,j+1  −→  vim+ j vim+ j+1 vim+ j+2 
Ii+1,j−1 Ii+1,j Ii+1,j+1 v ( i +1) m + j v ( i +1) m + j +1 v ( i +1) m + j +2

Obsérvese que las componentes de v I que aparecen para Iij se agru-

pan de 3 en 3 (en este caso, si el núcleo tiene tamaño 3 × 3): cada
grupo de 3 está separado una distancia m del grupo siguiente y el
grupo central corresponde a

(vim+ j , vim+ j+1 , vim+ j+2 ).

Por tanto, la componente k del vector imagen u J = Mv TI tiene la si-

guiente expresión, si k − 1 = lm + s

uk =K11 vk−m−1 + K12 vk−m + K13 vk−m+1

+K21 vk−1 + K22 vk + K23 vk+1
+K31 vk+m−1 + K32 vk+m + K33 vk+m+1

así que en la fila k de la matriz K todos los elementos son 0 excepto

los 9 que corresponden a los índices k − m − 1, k − m, k − m + 1, k −
1, k, k + 1 y k + m − 1, k + m, k + m + 1. Como es lógico, esto solo tiene
sentido para k > m y k < mn − m. El resto son 0.
Como estamos suponiendo una imagen de 5 × 4, la dimensión del
espacio inicial es 20. La matriz que se obtiene, está en la Figura 5. Se
explica a continuación. Téngase en cuenta que el ejemplo es de una
imagen muy pequeña (5 × 4), lo cual hace que todos los grupos de
entradas no nulas estén muy cerca. Como se ve, se ha multiplicado
todo por 1/40 para no poner decimales. Se observa lo siguiente:
• Todas las filas no nulas suman 1: esto es normal en trans-
formaciones digitales que “distribuyen la información” sin
reducirla.
• Las 5 primeras filas son todas de 0: esto se debe a que es-
tamos considerando el caso más fácil, en que el borde de la
imagen se descarta (lo que corresponde a la función nucleo1
de la Figura 2. Estas cinco filas son: los píxeles del borde su-
perior (hay 4) y el píxel del borde izquierdo de la fila 2 (estos
píxeles van seguidos en el vector v I ).
• Cada fila no nula contiene tres grupos de tres elementos se-
parados: estos grupos, 1, 2, 1 (que en realidad es 0.025, 0.05, 0.025),
 1. ALGUNOS EJEMPLOS IMPORTANTES 43

 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0 0
 
0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
1 
0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0
40 
0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0
 
0 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

FIGURA 5. Matrix 20 × 20 para el ejemplo 8.

2, 28, 2 y 1, 2, 1 son las filas de la matriz de la convolución K:

 
1 2 1
K = 1/40 2 28 2 .
1 2 1
• Cada dos píxeles con entradas no nulas, hay otros dos con
entradas nulas: corresponden a los bordes derecho e izquier-
do de cada fila y la siguiente.
• Como en la parte superior, en la inferior también hay 5 filas
llenas de 0.
Finalmente, si lo que hiciéramos fuera incluir los bordes exteriores
como “ceros”, la matriz que se obtendría sería la de la Figura 7 (se
deja como un ejercicio para el lector avezado). Esta transformación es
más útil que la anterior porque los bordes de la imagen no desapare-
cen sin más: se supone que la imagen inicial tiene un borde externo
de 0 y solo se pierde una pequeña parte de la información (como
se ve, ahora las filas correspondientes a puntos del borde no suman
1). Es importante observar cómo la diagonal principal tiene un valor
mayor que la suma del resto de los valores de la fila (y/o columna)
correspondiente (a este tipo de matrices se les llama dominantes por
44 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

−→

FIGURA 6. La imagen 5 × 4 convertida en una 3 × 2 con

borde nulo al usar un núcleo de tamaño 3 × 3 y omitir
los bordes.

 
28 2 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 
 
 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 
 
 0 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0 0 0 0 
 
 0 0 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0 0 0 
 
 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 0 0 
 
 0 
 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0 0 0 
 0 
 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0 0 
 0 
 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 0 
1 
 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0

0 
 
40  0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 0 
 
 0 0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 0 
 
 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 0 
 
 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 1 2 1 0 
 
 0 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 0 0 0 
 
 0 0 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 0 0 
 
 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 0 
 
 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 0 
 
 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 2 28 2 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 2 28

FIGURA 7. Matriz de cambio si se incluyen “bordes ex-

ternos nulos”.

filas/columnas). Para ejemplificarlo gráficamente, la Figura 8 mues-

tra cómo se convertiría una imagen 5 × 4 en escala de grises (64 to-
nos) utilizando el filtro de este ejemplo. Las imágenes en sus valores
 1. ALGUNOS EJEMPLOS IMPORTANTES 45

numéricos son:
   
35 36 41 36 31.25 38.65 41.75 33.27
39 42 38 33  30.92
 35.47 42.20 39.15 
45 45 44 34 −→ 40.07 44.92 43.20 31.60
   
50 47 47 30 44.37 46.82 45.52 28.82
53 51 48 33 46.15 48.92 42.02 28.42
y el núcleo que se utiliza, como ya se ha dicho, es
 
1 2 1
1 
2 28 2 .
40 1 2 1

−→

FIGURA 8. Conversión de una imagen 5 × 4 con el filtro

del ejemplo (un desenfoque incluyendo bordes nulos).

1.2. Análisis de tráfico. Un barrio de una ciudad puede presen-

tar, esquemáticamente, la estructura de calles que se muestra en la
Figura 9. Los números indican el flujo de coches por hora durante

70
250
x0,2 x1,2

y1,0 y1,1 y1,2

x0,1 x1,1
300 x31
y0,0 y0,1 y0,2
x0,0 x1,0
50
150

FIGURA 9. Tráfico esquemático de un barrio. Las varia-

bles xij son para el tráfico horizontal y las yij para el
vertical.
la hora punta del mediodía. Se quiere tener una idea del comporta-
miento del tráfico y de si se puede gestionar de alguna manera para
que en algún punto sea más fluido. Las variables xij indican, como
se explica en la figura, el tráfico horizontal por cada calle y las yij el
46 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

vertical. Se supone que el sentido de las calles es el que indican las

flechas (quizás poco realista), así que cada calle del modelo es de
sentido único (si hubiera calles de doble sentido, bastaría añadir un
eje en el grafo con una variable más).
Solo hay una ley que regula el flujo de tráfico: la de Kirchoff para
nodos: en cada nodo del grafo, el tráfico entrante ha de ser igual al
saliente. Por tanto, en el ejemplo, hay 9 requerimientos, cada uno de
los cuales impone una condición lineal (una suma igual a otra). En
concreto:
150 = x00 + y00
300 + y00 = x01 + y10
y10 = x02
x00 = x10 + y01
x01 + y01 = x11 + y11
x02 + y11 = x12 + 70
x10 = 50 + y02
x11 + y02 = x31 + y12
x12 + y12 = 250
La matriz aumentada asociada a este sistema en el orden de variables
indicado, es:
 
x00 x01 x02 x10 x11 x12 x31 y00 y01 y02 y10 y11 y12
 1 0 0 0 0 0 0 1 0 0 0 0 0 150 
 
 0 1 0 0 0 0 0 −1 0 0 1 0 0 300 
 
 0 0 1 0 0 0 0 0 0 0 −1 0 0 0
 
 −1 0 0 1 0 0 0 0 1 0 0 0 0 0
 
 0 −1 0 0 1 0 0 0 −1 0 0 1 0 0
 
 0 0 −1 0 0 1 0 0 0 0 0 −1 0 −70 
 0 0 0 −1 0 0 0 0 0 1 0 0 0 −50 
 
 0 0 0 0 −1 0 1 0 0 −1 0 0 1 0
0 0 0 0 0 −1 0 0 0 0 0 0 −1 −250
Incluso para un problema tan simple (piénsese lo complicado que
puede ser el tráfico de una ciudad entera), la dimensión del sistema
es de 9 ecuaciones y 13 variables. Este sistema resulta ser compatible
indeterminado pero incluso así, ya da algo de información. Por ejem-
plo (y este es un mero ejemplo de los muchos análisis que pueden
hacerse): las variables y12 , x00 y x11 pueden tomar cualquier valor. Por
tanto, en principio, cerrar esas tres calles podría seguir permitiendo el
flujo de coches que exigen las entradas y salidas impuestas. Sin em-
bargo, si se exige que las tres sean 0 (i.e. se cierran las tres calles al
tráfico), resulta que la solución del sistema da x00 = 0, y01 = −130,
lo cual significa que habría que cerrar la calle de x00 y cambiar de sentido
 1. ALGUNOS EJEMPLOS IMPORTANTES 47

la de y01 (esto, para un ayuntamiento, es peccata minuta, como todos

sabemos.
Se observa (como es natural) que la matriz del sistema es muy
dispersa: la mayor parte de las entradas son 0. Esto es consecuencia
de que las ecuaciones modelan un sistema real en que las relaciones
son locales (i.e. las condiciones en un nodo son independientes de las
de los nodos lejanos). Esta propiedad es bastante común en muchas
familias de problemas reales.
Se ha indicado ya que este ejemplo es un caso particular de las
Leyes de Kirchoff generales (hay una para nodos y otra para ciclos).
En el caso del flujo de tráfico, solo aplica la de los nodos. En circuitos
eléctricos aplican las dos.
1.3. Estática de un sistema de fuerzas: un puente. Supongamos
que se quiere construir un puente con la estructura de la Figura 10.
La distribución de fuerzas es la dada por las variables xi , yi , zi y los
y1 y2 y3

z0 z1 z2 z3 z4 z5 z6 z7
x0 x1 x2 x3

20 50 40

FIGURA 10. Esquema de un puente.

apoyos se sabe que están sujetos a las fuerzas indicadas (los dos ex-
tremos inferiores son fijos e inmóviles). Si se quiere que el sistema sea
estático, la suma de fuerzas en cada nodo ha de ser 0 (i.e. los nodos
no se mueven). Los ángulos se suponen de π/3 (60 grados). Se quiere
saber cuál es la distribución de esfuerzos (un esfuerzo negativo es
una compresión, uno positivo es una tensión) en los elementos.
Si el sistema es incompatible, el puente no se puede construir con
esas especificaciones. Si el sistema de ecuaciones es compatible inde-
terminado, el sistema (real) se denomina “hiperestático”. Si el sistema
de ecuaciones es compatible determinado, el sistema (real) se deno-
mina “isostático”.
Ha de tenerse en cuenta que, en cada nodo, hay que estudiar las
componentes horizontal y vertical de las fuerzas. Si los nodos supe-
riores son I1 , I2 , I3 , I4 y los inferiores (sin contar los extremos) J1 , J2 , J3 ,
entonces:
48 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

(1) La estabilidad del nodo I1 da

1 1
y1 − z0 + z1 = 0 (horizontal)
√ 2 √2
3 3
z0 + z1 = 0 (vertical)
2 2
Como hay 7 nodos “efectivos” (los extremos, al estar sujetos al te-
rreno, no añaden ni quitan esfuerzos), en total se obtienen 14 ecua-
ciones (que no vamos a escribir). Igual que en las de I1 , hay una con-
dición vertical y una horizontal en cada una, y cada ecuación (en la
estructura de la Figura 10) contiene a lo sumo 3 variables, por lo que
el sistema de ecuaciones que se obtiene es también muy disperso (ca-
da fila contiene 14 entradas, de la cuales por lo menos 11 son nulas).
Claramente, cuantos más elementos posee una estructura, mayor
es el número de variables y, usualmente, de ecuaciones, con lo que es
fácil obtener sistemas de decenas de ecuaciones e incógnitas.
Si, además, hubiera más fuerzas actuando fuera de los nodos, ha-
bría que utilizar también el principio de “anulación de momentos”
(que no aplica en este ejemplo).

2. El algoritmo de Gauss y la factorizacón LU

Partiendo de la Ecuación (1), el algoritmo de Gauss consiste en
transformar A y b mediante operaciones “sencillas” en una matriz Ã
triangular superior y un vector b̃, para que el nuevo sistema Ãx = b̃
sea directamente soluble mediante sustitución regresiva (es decir, se
puede calcular la variable xn despejando directamente y, sustituyen-
do “hacia arriba” en la ecuación n − 1, se puede calcular xn−1 , etc.) Es
obvio que se requiere que la solución del sistema nuevo sea la misma
que la del original. Para ello, se permite realizar la siguiente opera-
ción:
• Se puede sustituir una ecuación Ei (la fila i −ésima de A) por
una combinación lineal de la forma Ei + λEk donde k < i y
λ ∈ R. Si se hace esto, también se sustituye bi por bi + λbk .
El hecho de que la combinación tenga coeficiente 1 en Ei es lo que
obliga a que las soluciones del sistema modificado sean las mismas
que las del original.
LEMA 1. Para transformar una matriz A en una matriz à según la ope-
ración anterior, basta multiplicar A por la izquierda por la matriz Lik (λ)
cuyos elementos son:
• Si m = n, entonces ( Lik (λ))mn = 1 (diagonal de 1).
 2. EL ALGORITMO DE GAUSS Y LA FACTORIZACÓN LU 49

• Si m = i, n = k, entonces ( Lik (λ))mn = λ (el elemento (i, k) es

λ).
• Cualquier otro elemento es 0.
EJEMPLO 10. Si se parte de la matriz A
 
3 2 −1 4
0 1 4 2 
A= 6 −1 2

5 
1 4 3 −2
y se combina la fila 3 con la fila 1 multiplicada por −2 (para “hacer
un cero” en el 6), entonces se ha de multiplicar A por L31 (−2)
    
1 0 0 0 3 2 −1 4 3 2 −1 4
 0 1 0 0 0 1 4 2   2 
   = 0 1 4 .
 −2 0 1 0 6 −1 2 5  0 −5 4 −3
0 0 0 1 1 4 3 −2 1 4 3 −2
De manera simplificada, el algoritmo de Gauss puede enunciarse
como indica el Algoritmo 7.
La línea marcada con el asterisco en el Algoritmo 7 es exactamente
la multiplicación de à por la izquierda por la matriz L ji (m ji ), donde
el multiplicador m ji es justamente
m ji = − A ji /Aii .
Así que al final, la matriz Ã, triangular superior, es un producto de
esas matrices:
(2) Ã = Ln,n−1 (mn,n−1 ) Ln,n−2 (mn,n−2 ) · · · L2,1 (m2,1 ) A = L̃A
donde L̃ es una matriz triangular inferior con 1 en la diagonal (esto
es un sencillo ejercicio). Resulta también sencillo comprobar que (y
esto realmente parece magia) que
LEMA. La matriz inversa del producto de todas las matrices de (2) es la
matriz triangular inferior que en cada componente ( j, i ) contiene el valor
−m ji .
Así que, sin más complicación, se ha demostrado el siguiente re-
sultado:
TEOREMA 3. Si en el proceso de reducción de Gauss no hay ningún ele-
mento de la diagonal igual a 0, entonces existe una matriz L triangular infe-
rior cuyos elementos son los sucesivos multiplicadores cambiados de signo
en su posición correspondiente y una matriz diagonal superior U tal que
A = LU
50 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

Algoritmo 7 (Algoritmo de Gauss para sistemas lineales)

Input: Una matriz A y un vector b, ambos de orden n
Output: O bien un mensaje de error o bien una matriz à y un vector
b̃ tales que à es triangular superior y que el sistema Ãx = b̃ tiene
las mismas soluciones que el Ax = b
⋆INICIO
à ← A, b̃ ← b, i ← 1
while i < n do
if Ãii = 0 then
return ERROR [división por cero]
end if
[combinar cada fila bajo la i con la i]
j ← i+1
while j ≤ n do
m ji ← − Ã ji / Ãii
[La siguiente linea es un bucle, cuidado]
à j ← à j + m ji Ãi [*]
b̃ j ← b̃ j + m ji b̃i
j ← j+1
end while
i ← i+1
end while
return Ã, b̃

y tal que el sistema Ux = b̃ = L−1 b es equivalente al sistema incial Ax = b.

Con este resultado, se obtiene una factorización de A que simplifi-

ca la resolución del sistema original, pues se puede reescribir Ax = b
como LUx = b; yendo por partes, se hace:

• Primero se resuelve el sistema Ly = b, por sustitución directa

—es decir, de arriba abajo, sin siquiera dividir.
• Luego se resuelve el sistema Ux = y, por sustitución regresiva
—es decir, de abajo arriba.

Esta manera de resolver solo requiere conservar en memoria las ma-

trices L y U y es muy rápida por comparación con el método de Cra-
mer, por ejemplo, o como se explica a continuación, con el uso de la
inversa.
 2. EL ALGORITMO DE GAUSS Y LA FACTORIZACÓN LU 51

2.1. [Contenido no esencial] Complejidad del algoritmo de Gauss.

Clásicamente, la complejidad de los algoritmos que utilizan opera-
ciones aritméticas se medía calculando el número de multiplicacio-
nes necesarias (pues la multiplicación era una operación mucho más
compleja que la adición). Hoy día esta medida es menos representati-
va, pues las multiplicaciones en coma flotante se realizan en un tiem-
po prácticamente equivalente al de las adiciones (con ciertos matices,
pero este es uno de los avances importantes en velocidad de proce-
samiento).
De esta manera, si se supone que las divisiones pueden hacerse con
exactitud1, se tiene el siguiente resultado:
LEMA 2. Si en el proceso de simplificación de Gauss no aparece ningún
cero en la diagonal, tras una cantidad de a lo sumo (n − 1) + (n − 2) +
· · · + 1 combinaciones de filas, se obtiene un sistema Ãx = b̃ donde à es
triangular superior.
Cada combinación de filas en el paso k, tal como se hace el algo-
ritmo requiere una multiplicación por cada elemento de la fila (des-
contando el que está justo debajo del pivote, que se sabe que hay que
sustituir por 0), que en este paso da n − k productos, más una divi-
sión (para calcular el multiplicador de la fila) y otra multiplicación
para hacer la combinación del vector. Es decir, en el paso k hacen falta
( n − k )2 + ( n − k ) + ( n − k )
operaciones “complejas” (multiplicaciones, esencialmente) y hay que
sumar desde k = 1 hasta k = n − 1, es decir, hay que sumar
n −1 n −1
∑ i2 + 2 ∑ i.
i =1 i =1

La suma de los primeros r cuadrados es (clásica) r (r + 1)(2r + 1)/6,

mientras que la suma de los primeros r números es r (r + 1)/2. Así
pues, en total, se tienen
( n − 1) n (2( n − 1) + 1) n3 n2 5n
+ ( n − 1) n = + − .
6 3 2 6
Para ahora resolver el sistema triangular superior Ãx b̃, hace falta una
división por cada línea y k − 1 multiplicaciones en la fila n − k (con
k desde 0 hasta n), así que hay que sumar n + 1 + · · · + (n − 1) =

1Lo cual es, como ya se ha dicho muchas veces, demasiado suponer.

52 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

n ( n +1)
2 .
Por tanto, para resolver un sistema con el método de Gauss,
hacen falta en total
n3 n
+ n2 − operaciones para Ax = b.
3 3
Si el algoritmo se utiliza para resolver m sistemas de la forma
Ax = bi (para diferentes bi ), todas las operaciones son iguales para
triangular A y simplemente hay que recomputar b̃i y resolver, Esto
requiere (n − 1) + · · · + 1 = n(n − 1)/2 multiplicaciones (las del pro-
ceso de triangulación) y 1 + 2 + · · · + n para “despejar”. Descontan-
do las que se hicieron en la resolución de b, resulta que para resolver
los m sistemas, hacen falta:
n3 n
(3) + mn2 − operaciones para m sistemas.
3 3
2.1.1. Comparación con utilizar A−1 . Es sencillo comprobar que el
cálculo general de la inversa de una matriz requiere, utilizando el
algoritmo de Gauss-Jordan (o por ejemplo, resolviendo los sitemas
Ax = ei para la base estándar) al menos n3 operaciones. Hay mejores
algoritmos (pero se están intentando comparar métodos análogos).
Una vez calculada la inversa, resolver un sistema Ax = b consiste en
multiplicar b a la izquierda por A−1 , que requiere (obviamente) n2
productos. Por tanto, la resolución de m sistemas requiere
n3 + mn2 operaciones complejas.
Que siempre es más grande que (3). Así que, si se utiliza el método
de Gauss, es mejor conservar la factorización LU y utilizarla para “la
sustitución” que calcular la inversa y utilizarla para multiplicar por
ella.
Claro está que esta comparación es entre métodos análogos: hay mane-
ras de computar la inversa de una matriz en menos (bastante menos)
de n3 operaciones (aunque siempre más que un orden de n2 log(n)).
2.2. Estrategias de Pivotaje, el algoritmo LUP. Como ya se dijo,
si en el proceso de reducción gaussiana aparece un pivote (el elemen-
to que determina el multiplicador) con valor 0 (o incluso con deno-
minador pequeño), o bien puede no continuarse de manera normal o
bien puede que aparezcan errores muy grandes debido al redondeo.
Esto puede aliviarse utilizando estrategias de pivotaje: cambiando el
orden de las filas o de las columnas. Si solo se cambia el orden de las
filas, se dice que se realiza un pivotaje parcial. Si se hacen ambas ope-
raciones, se dice que se realiza un pivotaje total. En estas notas solo se
estudiará el pivotaje parcial.
 2. EL ALGORITMO DE GAUSS Y LA FACTORIZACÓN LU 53

DEFINICIÓN 9. Una matriz de permutación es una matriz cuadrada

formada por ceros salvo que en cada fila y en cada columna hay exac-
tamente un 1.
(Para que una matriz sea de permutación basta con que se cons-
truya a partir de la matriz identidad, permutando filas).
Es obvio que el determinante de una matriz de permutación es
distinto de 0 (de hecho, es bien 1 bien −1). No es tan sencillo com-
probar que la inversa de una matriz de permutación P es también
una matriz de permutación y, de hecho, es P T (su traspuesta).
LEMA 3. Si A es una matriz n × m y P es una matriz cuadrada de
orden n, de permutación, que solo tiene 2 elementos fuera de la diagonal
no nulos, digamos los (i, j) y ( j, i ) (pues tiene que ser simétrica), entonces
PA es la matriz obtenida a partir de A intercambiando las filas i y j. Para
intercambiar por columnas se ha de hacer a la derecha (pero no se explica
aquí).
DEFINICIÓN 10. Se dice que en el proceso de reducción de Gauss
se sigue una estrategia de pivotaje parcial si el pivote en el paso i es
siempre el elemento de la columna i de mayor valor absoluto.
Para realizar el algoritmo de Gauss con pivotaje parcial, basta rea-
lizarlo paso a paso salvo que, antes de fijar el pivote, se ha de buscar,
bajo la fila i el elemento de mayor valor absoluto de la columna i. Si
este está en la fila j (y será siempre j ≥ i, al buscar por debajo de i),
entonces se han de intercambiar las filas i y j de la matriz.
La estrategia de pivotaje da lugar a una factorización que no es
exactamente la LU, sino con una matriz de permutación añadida:
LEMA 4. Dado el sistema de ecuaciones Ax = b con A no singular,
siempre existe una matriz P de permutación y dos matrices L, U, la primera
triangular inferior, la segunda triangular superior, con
PA = LU.
La prueba de este resultado no es directa, hay que hacer un razo-
namiento por recuerrencia (sencillo pero no inmediato).
Tanto Octave como Matlab incluyen la función lu que, dada una
matriz A, devuelve tres valores: L, U y P.
Por ejemplo, si
 
1 2 3 4
 −1 −2 5 6
A=  −1 −2 −3 7

0 12 7 8
54 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

entonces
     
1 0 0 0 1 2 3 4 1 0 0 0
 0 1 0 0  
0 12 7 8   1
L= U= 0 0 0 
 −1 0 1 0 0 0 8 10 , P = 0 1 0 0
−1 0 0 1 0 0 0 11 0 0 1 0
Para calcular L, U y P no hay más que realizar el algoritmo ordi-
nario de Gauss salvo que cuando se realice un intercambio entre la fila i
y la fila j se ha de realizar el mismo en la L hasta entonces calculada
(cuidado: solo en las primeras i − 1 columnas, la zona de la diagonal
con “1” no se ha de tocar) y multiplicar la P ya calculada (comenzan-
do con P = Idn ) por la izquierda por Pij (la matriz de permutación
de las filas i y j).
En fin, puede expresarse el algoritmo LUP como en el Algoritmo
8.
2.3. El número de condición: comportamiento del error relati-
vo. ¿Cómo es de “estable” la resolución de un sistema de ecuaciones
lineales? Una manera muy simple de acercarse a este problema es com-
parar los “errores relativos” en una solución si se cambia el vector b
por uno modificado un poco. Supóngase que en lugar del sistema (1)
original, cuya solución es x, se tiene uno modificado
Ay = b + δb
donde δb es un vector pequeño. La solución será de la forma x + δx,
para cierto δx (que uno espera que sea pequeño).
La manera de medir tamaños de vectores es mediante una norma
(la más común es la longitud, en el espacio eculídeo, pero no es la que
se utilizará aquí). Como x es una solución, se tiene que
A( x + δx ) = b + δb,
así que
Aδx = δb,
pero como se parte de una modificación de b y estudiar cómo se mo-
difica x, despejando y queda
δx = A−1 δb
y midiendo tamaños (es decir, normas, que se denotan ∥∥), quedaría
∥δx ∥ = ∥ A−1 δb∥.
Se está estudiando el desplazamiento relativo, que es más relevante que
el absoluto. Ahora se ha incluir ∥ x ∥ en el primer miembro. La infor-
mación que se tiene es que Ax = b, de donde ∥ Ax ∥ = ∥b∥. Así que
 2. EL ALGORITMO DE GAUSS Y LA FACTORIZACÓN LU 55

Algoritmo 8 (Factorización LUP para una matriz A)

Input: Una matriz A de orden n
Output: O bien un mensaje de error o bien tres matrices: L triangu-
lar inferior con 1 en la diagonal, U triangular superior y P matriz
de permutación tales que LU = PA
⋆INICIO
L ← Idn , U ← A, P ← Idn , i ← 1
while i < n do
p ← fila tal que |U pi | es máximo, con p ≥ i
if U pi = 0 then
return ERROR [división por cero]
end if
[intercambiar filas i y p]
P ← Pip P
U ← Pip U
[en L solo se intercambian las filas i y p de la submatriz
n × (i − 1) de la izquierda, ver texto]
L ← L̃
[combinar filas en U y llevar cuenta en L]
j ← i+1
while j <= n do
m ji ← Uji /Uii
Uj ← Uj − mij Ui
L ji ← m ji
j ← j+1
end while
i ← i+1
end while
return L, U, P

queda

∥δx ∥ ∥ A−1 δb∥

= ,
∥ Ax ∥ ∥b∥

pero esto no dice mucho (pues es una identidad obvia, se querría

acotar el desplazamiento relativo de la solución si se desplaza un poco
el término independiente).
56 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

Supóngase que existe determinado objeto llamado norma de una

matriz, que se denotaría ∥ A∥ que cumple que ∥ Ax ∥ ≤ ∥ A∥∥ x ∥. En-
tonces el miembro de la izquierda la ecuación anterior, leída de de-
recha a izquierda quedaría
∥ A−1 δb∥ ∥δx ∥ ∥δx ∥
= ≥
∥b∥ ∥ Ax ∥ ∥ A∥∥ x ∥
y, por otro lado, aplicando el mismo razonamiento a la parte “de las
b”, se tendría que
∥ A−1 ∥∥δb∥ ∥ A−1 δb∥
≥ ,
∥b∥ ∥b∥
y combinando todo,
∥ A−1 ∥∥δb∥ ∥δx ∥
≥
∥b∥ ∥ A∥∥ x ∥
de donde, finalmente, se obtiene una cota superior para el desplaza-
miento relativo de x:
∥δx ∥ ∥δb∥
≤ ∥ A∥∥ A−1 ∥
∥x∥ ∥b∥
El caso es que tal objeto, llamado norma de una matriz (o más bien de
una aplicación lineal) existe. De hecho, existen muchos. En estas notas
se va a utilizar el siguiente:
DEFINICIÓN 11. Se denomina norma infinito de una matriz cuadra-
da A = ( aij ) al número
n
∥ A∥∞ = max ∑ |aij |,
1≤ i ≤ n j =1

es decir, el máximo de las sumas de los valores absolutos por filas.

En realidad, se podrían utilizar muchas otras definiciones, pero
esta es la que se usará en estas notas.
LEMA 5. La norma infinito cumple que, para cualquier vector x, ∥ Ax ∥∞ ≤
∥ A∥∞ ∥ x ∥∞ , donde ∥ x ∥∞ es la norma dada por el máximo de los valores ab-
solutos de las coordenadas de x.
Es decir, si se toma como medida del tamaño de un vector el máximo
de sus coordenadas en valor absoluto, y lo denominamos ∥ x ∥∞ se
tiene que
∥δx ∥∞ ∥δb∥∞
≤ ∥ A ∥ ∞ ∥ A −1 ∥ ∞ .
∥ x ∥∞ ∥b∥∞
 2. EL ALGORITMO DE GAUSS Y LA FACTORIZACÓN LU 57

Al producto ∥ A∥∞ ∥ A−1 ∥∞ se le denomina número de condición de la

matriz A para la norma infinito, se denota κ ( A) y es una medida del má-
ximo desplazamiento posible de la solución si se desplaza un poco el
vector. Así que, cuanto más grande sea el número de condición peor
se comportará en principio el sistema respecto de pequeños cambios
en el término independiente.
El número de condición también sirve para acotar inferiormente el
error relativo cometido:

LEMA 6. Sea A una matriz no singular de orden n y x una solución del

sistema Ax = b. Sea δb un “desplazamiento” de las condiciones iniciales y
δx el “desplazamiento” correspondiente en la solución. Entonces:

1 ∥δb∥ ∥δx ∥ ∥δb∥

≤ ≤ κ ( A) .
κ ( A) ∥b∥ ∥x∥ ∥b∥

Así que se puede acotar el error relativo con el residuo relativo (el número
∥δb∥/∥b∥).

EJEMPLO 11. Considérese el sistema

0.853x + 0.667y = 0.169

0.333x + 0.266y = 0.067

Su número de condición para la norma infinito es 376.59, así que un

cambio de una milésima relativo en las condiciones iniciales (el vec-
tor b) puede dar lugar a un cambio relativo de más del 37% en la
solución. La de dicho sistema es x = 0.055+, y = 0.182+. Pero el
número de condición tan grande avisa de que una pequeña pertur-
bación originará grandes modificaciones en la solución. Si se pone,
en lugar de b = (0.169, 0.067) el vector b = (0.167, 0.067) (un des-
plazamiento relativo del 1.1% en la primera coordenada), la solución
es x = −0.0557+, y = 0.321+. Por un lado, la x no tiene siquiera
el mismo signo, por otro, el desplazamiento relativo en la y es del
76%, totalmente inaceptable. Imagínese que el problema describe un
sistema estático de fuerzas y que las medidas se han realizado con
instrumentos de precisión de ±10−3 . Para b = (0.168, 0.067), la di-
ferencia es menor (aunque aun inaceptable en y), pero el cambio de
signo persiste…
58 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

2.4. Un ejemplo explícito de Gauss con pivotaje. Se muestra a

continuación el desarrollo completo del método de Gauss con pivo-
taje parcial aplicado a un sistema (matriz de coeficientes) 4 × 4. Sea
   
1 6 8 3 2
 2 4 6 0  
A= , b = 1
 −1 3 6 1 2
0 4 2 8 3

Se quieren calcular las matrices L (triangular inferior), U (triangular

superior) y P (matriz de permutación) que satisfagan:

PA = LU

pero vamos a llevar cuenta de b, también, así que en todo el desa-

rrollo trabajaremos con la matriz ampliada. Comenzamos, por tanto,
haciendo
 
1 6 8 3 2
 2 4 6 0 1
U=  −1 3 6 1 2 ,

0 4 2 8 3

que es la matriz que vamos a ir transformando hasta que la parte

de los coeficientes sea triangular superior. Procedemos a realizar el
algoritmo, paso a paso.
(1) Se han de intercambiar las filas 1 y 2 para que el pivote tenga
el valor 2. Para esto, se usa la matriz P12 , y se obtiene
    
2 4 6 0 1 0 1 0 0 1 6 8 3 2
 1 6 8 3 2 1 0 0 0  2 4 6 0 1
U1 ← 
 −1
=  ,
3 6 1 2 0 0 1 0  −1 3 6 1 2
0 4 2 8 3 0 0 0 1 0 4 2 8 3

(2) Utilizando las matrices L12 (−1/2) y L13 (1/2) se eliminan

los términos 1 y −1, respectivamente, de debajo del pivote
de la primera columna. Recuérdese que
   
1 0 0 0 1 0 0 0
− 1 1 0 0  0 1 0 0
L12 (−1/2) =  2
 0 0 , L (1/2) =  ,
1 0 13 1
2 0 1 0
0 0 0 1 0 0 0 1
 2. EL ALGORITMO DE GAUSS Y LA FACTORIZACÓN LU 59

y, haciendo todos los cálculos, queda

 
2 4 6 0 1
0 4 5 3 23 
U2 ← L13 (1/2) · L12 (−1/2)U1 =  0 5 9 1 5 

2
0 4 2 8 3
(3) Ahora hay que intercambiar las filas 2 y 3 para que el 5 quede
como pivote. Para esto se utiliza P23 :
 
2 4 6 0 1
0 5 9 1 25 
U3 ← P23 U2 =  0 4 5 3 3  .

2
0 4 2 8 3
(4) Se eliminan los dos 4 que hay debajo del pivote, utilizando
L23 (−4/5) y L24 (−4/5):
 
2 4 6 0 1
0 5 9 1 5
 2
U4 ←   = L24 (−4/5) L23 (−4/5)U3 .
0 0 − 11 5
11
5 − 12 
0 0 − 265
36
5 1
(5) Se intercambian las filas 3 y 4 para que el −26/5 quede como
pivote, utilizando P34 :
 
2 4 6 0 1
0 5 9 1 5
 2
U5 ←  .
0 0 − 26 5
36
5 1
0 0 − 11 5
11
5 − 12
(6) Y, finalmente, se elimina el término −11/5 utilizando L34 (−11/26):
 
2 4 6 0 1
 5
0 5 9 1 
U6 ←  0 0 − 26 36
2 .

 5 5 1 
0 0 0 − 13 − 13
11 12

Así pues, se puede transformar el sistema original en

 
2 4 6 0 1
0 5 9 1 5
 2
U6 =  ,
0 0 − 26 5
36
5 1 
0 0 0 − 1113 − 13
12
60 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

que es triangular superior, con la siguiente secuencia de matrices:

U6 = L34 (−11/26) P34 L24 (−4/5) L23 (−4/5) P23 L13 (1/2) L12 (−1/2) P12 A.

Para construir L, hay que utilizar los multiplicadores cambiados de

signo y hacer los intercambios parciales de filas según indiquen los
Pij , sin contar el primero:
(1) Se construye la matriz de multiplicadores con el signo cam-
biado para la primera columna:
 
1 0 0 0
 1 1 0 0
 
L1 =  21 
 − 2 0 1 0
0 0 0 1

(2) Se intercambia la parte anterior a la columna 2 según indique

el segundo cambio de filas, P23 (es decir, filas 2 y 3):
 
1 0 0 0
 − 1 1 0 0
 
L2 =  21 
 2 0 1 0
0 0 0 1

como se ve, la diagonal de 1 sigue estando ahí.

(3) Se introducen los multiplicadores cambiados de signo en la
segunda columna:
 
1 0 0 0
 1 
 − 2 1 0 0

L3 =  1 4 .
1 0 
 2 5 
4
0 5 0 1

(4) Se intercambia la parte izquierda de las filas correspondien-

tes a la siguiente matriz P, que es P34 , es decir, las filas 3 y
4:
 
1 0 0 0
 1 
 − 2 1 0 0
L4 = 
 0 4 1 0 .

 5 
1 4
2 5 0 1
 3. ALGORITMOS APROXIMADOS (DE PUNTO FIJO) 61

(5) Se termina introduciendo el último multiplicador, cambiado

de signo:
 
1 0 0 0
 1 
 − 2 1 0 0

L4 =  .
0 4
1 0 
 5 
1 4 11
2 5 26 1
Con esto, se obtiene que
P34 P23 P12 A = L4 U6 ,
donde es fácil ver que
 
0 1 0 0
0 0 1 0
P = P34 P23 P12 = 
0
.
0 0 1
1 0 0 0

3. Algoritmos aproximados (de punto fijo)

Tanto la complejidad computacional del algoritmo de Gauss co-
mo su inestabilidad respecto a la división por pequeños números
(debida a la introducción inevitable de errores de redondeo) llevan a
plantearse buscar otro tipo de algoritmos con mejores propiedades.
En algunas situaciones sería incluso utópico plantearse resolver un
sistema en un cantidad similar a n3 multiplicaciones (pongamos por
caso que n ≃ 106 , entonces n3 ≃ 1018 y haría falta demasiado tiempo
para realizar todas las operaciones).
Si se parte de un sistema como el original (1), de la forma Ax = b,
puede tratar de convertirse en un problema de punto fijo mediante
una separación de la matriz A en dos, A = N − P, donde N es in-
vertible. De este modo, el sistema Ax = b queda ( N − P) x = b, es
decir
( N − P) x = b ⇒ Nx = b + Px ⇒ x = N −1 b + N −1 Px,
si se llama c = N −1 b y M = N −1 P, queda el siguiente problema de
punto fijo
x = Mx + c
que, si puede resolverse, puede hacerse mediante iteraciones de la
misma forma que en el Capítulo 2: se comienza con una semilla x0 y
se itera
xn = Mxn−1 + c,
62 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

hasta alcanzar una precisión suficiente. Hacen falta los siguientes re-
sultados.

TEOREMA 4. Supongamos que M es una matriz de orden n y que ∥ M∥∞ <

1. Entonces la ecuación x = Mx + c tiene solución única para todo c y la
iteración xn = Mxn−1 + c converge a ella para cualquier vector inicial x0 .

TEOREMA 5. Dada la matriz M con ∥ M∥∞ < 1, y dado x0 una semilla

para el método iterativo del Teorema 4, si s es la solución del problema x =
Mx + c, se tiene la siguiente cota:

∥ M∥n∞
∥ xn − s∥∞ ≤ ∥ x1 − x0 ∥ ∞ .
1 − ∥ M ∥∞

Recuérdese que, para vectores, la norma infinito ∥ x ∥∞ viene dada por el

mayor valor absoluto de las componentes de x.

Con estos resultados, se explican los dos métodos básicos itera-

tivos para resolver sistemas de ecuaciones lineales: el de Jacobi, que
corresponde a tomar N como la diagonal de A y el de Gauss-Seidel,
que corresponde a tomar N como la parte triangular inferior de A
incluyendo la diagonal.

3.1. El algortimo de Jacobi. Si en el sistema Ax = b se “despeja”

cada coordeanada xi , en función de las otras, queda una expresión así:

1
xi = (b − ai1 x1 − · · · − aii−1 xi−1 − aii+1 xi+1 − · · · − ain xn ),
aii i

en forma matricial,

x = D −1 ( b − ( A − D ) x ),

donde D es la matriz diagonal cuyos elementos no nulos son exacta-

mente los de la diagonal de A. En fin, puede escribirse, por tanto,

x = D −1 b − D −1 ( A − D ) x,

que es una ecuación de tipo punto fijo. Si la matriz D −1 ( A − D ) cum-

ple las condiciones del Teorema 4, entonces la iteración de Jacobi con-
verge para cualquier semilla x0 y se tiene la cota del Teorema 5. Para
comprobar las condiciones, hace falta calcular D −1 ( A − D ), aunque
puede comprobarse de otras maneras (véase el Lema 7).
 4. ANEXO: CÓDIGO EN MATLAB/OCTAVE 63

3.2. El algoritmo de Gauss-Seidel. Si en lugar de utilizar la dia-

gonal de la matriz A para descomponerla, se utiliza la parte trian-
gular inferior (incluyendo la diagonal), se obtiene un sistema de la
forma
x = T −1 b − T −1 ( A − T ) x,
que también es una ecuación de punto fijo. Si la matriz T −1 ( A − T )
cumple las condiciones del Teorema 4, entonces la iteración de Gauss-
Seidel converge para cualquier semilla x0 , y se tiene la cota del Teo-
rema 5. Para comprobar las condiciones, haría falta calcular T −1 ( A −
T ), aunque puede comprobarse de otras maneras.
DEFINICIÓN 12. Se dice que una matriz A = ( aij ) es estrictamente
diagonal dominante por filas si
| aii | > ∑ |aij |
j ̸ =i

para todo i entre 1 y n.

Para estas matrices, tanto el método de Jacobi como el de Gauss-
Seidel son convergentes:
LEMA 7. Si A es una matriz estrictamente diagonal dominante por
filas, entonces tanto el método de Jacobi como el de Gauss-Seidel convergen
para cualquier sistema de la forma Ax = b.
Para el de Gauss-Seidel, además, se tiene que
LEMA 8. Si A es una matriz simétrica definida positiva, entonces el mé-
todo de Gauss-Seidel converge para cualquier sistema de la forma Ax = b.

4. Anexo: Código en Matlab/Octave

Se incluye el código de varios de los algoritmos descritos, usable
tanto en Matlab como en Octave.
4.1. El algoritmo de Gauss. El siguiente código implementa el
algoritmo de reducción de Gauss para un sistema Ax = b y devuelve
la matriz L, la matriz U y el vector b transformado, suponiendo que los
multiplicadores por defecto nunca son 0. Si alguno es cero, se abandona
el procedimiento.
La entrada ha de ser:
A: una matriz cuadrada (si no es cuadrada, la salida es la trian-
gulación por la diagonal principal),
b: un vector con tantas filas como columnas tiene A.
La salida es una terna L, At, bt, como sigue:
 64 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

function [L, At , bt] = gauss(A,b)

n = size(A);
m = size(b);
4 if(n(2) ~= m(1))
warning('The sizes of A and b do not match ');
return;
end
At=A; bt=b; L=eye(n);
k=1;
while (k<n(1))
l=k+1;
if(At(k,k) == 0)
warning('There is a 0 on the diagonal ');
14 return;
end
% careful with rows & columns:
% L(l,k) means ROW l, COLUMN k
while(l<=n)
L(l,k)=At(l,k)/At(k,k);
% Combining rows is easy in Matlab
At(l,k:n) = [0 At(l,k+1:n) - L(l,k) * At(k,k+1:n)];
bt(l)=bt(l)-bt(k)*L(l,k);
l=l+1;
24 end
k=k+1;
end
end

FIGURA 11. Código del algoritmo de reducción de Gauss

L: es la matriz triangular inferior (de los multiplicadores),

At: es la matriz A transformada (la reducida), que se denomi-
na U en la factorización LU, y es triangular superior.
bt: el vector transformado.
De manera, que el sistema que se ha de resolver, equivalente al ante-
rior, es At × x = bt.
4.2. El algoritmo LUP. Como se vio arriba, el algoritmo de re-
ducción de Gauss está sujeto a que no aparezca ningún cero como
pivote y además puede dar lugar a errores importantes de redondeo
si el pivote es pequeño. A continuación se muestra un código en Ma-
talb/Octave que implementa el algoritmo LUP, por el cual se obtiene
una factorización LU = PA, donde L y U son triangular inferior y
superior, respectivamente, la diagonal de L está formada por 1 y P es
una matriz de permutación. La entrada es:
A: una matriz cuadrada de orden n.
b: un vector de n filas.
 4. ANEXO: CÓDIGO EN MATLAB/OCTAVE 65

Devuelve cuatro matrices, L, At, P y bt, que corresponded a L, U,

P y al vector transformado según el algoritmo.
 66 3. SOLUCIÓN APROXIMADA DE SISTEMAS LINEALES

function [L, At , P, bt] = gauss_pivotaje(A,b)

n = size(A);
3 m = size(b);
if(n(2) ~= m(1))
warning('Dimensions of A and b do not match ');
return;
end
At=A;
bt=b;
L=eye(n);
P=eye(n);
i=1;
13 while (i<n)
j=i+1;
% beware nomenclature:
% L(j,i) is ROW j, COLUMN i
% the pivot with greatest absolute value is sought
p = abs(At(i,i));
pos = i;
for c=j:n
u = abs(At(c,i));
if(u>p)
23 pos = c;
p = u;
end
end
if(u == 0)
warning('Singular system ');
return;
end
% Swap rows i and p if i != p
% in A and swap left part of L
33 % This is quite easy in Matlab , there is no need
% for temporal storage
P([i pos],:) = P([pos i], :);
if(i ~= pos)
At([i pos], :) = At([pos i], :);
L([i pos], 1:i-1) = L([pos i], 1:i-1);
b([i pos], :) = b([pos i], :);
end
while(j<=n)
L(j,i)=At(j,i)/At(i,i);
43 % Combining these rows is easy
% They are 0 up to column i
% And combining rows is easy as above
At(j,i:n) = [0 At(j,i+1:n) - L(j,i)*At(i,i+1:n)];
bt(j)=bt(j)-bt(i)*L(j,i);
j=j+1;
end
i=i+1;
end
end

FIGURA 12. Código del algoritmo LUP

 CAPÍTULO 4

Interpolación

Dada una colección de datos, la tentación humana es utilizarlos

como fuentes de conocimiento en lugares desconocidos. Específica-
mente, dada una lista de coordenadas ligadas a un tipo de suceso
(un experimento, unas mediciones…) ( xi , yi ), lo “natural” es intentar
utilizarla para deducir o predecir el valor que tomaría la y si la x fuera
cualquier otra. Este es el afán interpolador y extrapolador del hombre.
No se puede hacer nada para evitarlo. Lo que se puede hacer es calcu-
lar las maneras más razonables de llevar a cabo dicha interpolación.
En todo este tema se parte de una lista de datos
x x0 x1 . . . x n −1 x n
(4)
y y0 y1 . . . y n −1 y n
que se supone ordenada en las coordenadas x, que además son dife-
rentes: xi < xi+1 . El objetivo es encontrar funciones que —de alguna
manera— tengan relación (cierta cercanía) con dicha lista de datos o
nube de puntos.

1. Interpolación lineal (a trozos)

La primera idea (sencilla pero funcional) es utilizar la función
definida a trozos entre x0 y xn que consiste en “unir cada punto con el
siguiente con una recta”. Esto se denomina interpolación lineal a trozos
ó spline lineal (se definirá spline más adelante con generalidad).
DEFINICIÓN 13. La función de interpolación lineal a trozos de la tabla
(4) es la función f : [ x0 , xn ] → R definida de la siguiente manera:
y − y i −1
f (x) = i ( x − xi−1 ) + yi−1 si x ∈ [ xi−1 , xi ]
x i − x i −1
es decir, la función definida a trozos que consiste en los segmentos
que unen ( xi−1 , yi−1 ) con ( xi , yi ), definida desde x0 hasta xn , para
i = 1, . . . , n.
La interpolación lineal tiene varias características que la hacen in-
teresante:
• Es muy sencilla de calcular.
67
68 4. INTERPOLACIÓN

f (x)
interpolación
4

4 6 8 10

FIGURA 1. Interpolación lineal de una nube de 9 puntos.

Compárese con la función original (en verde).

• Pasa por todos los puntos de la tabla de datos.

• Es continua.
Por eso se utiliza con frecuencia para representar funciones (es lo
que hace Matlab por defecto), pues si la nube de puntos es densa, los
segmentos serán pequeños y las esquinas se notarán poco.
La pega de la interpolación lineal son precisamente las esquinas
que aparecen siempre que la tabla de datos no corresponda a puntos
de una recta.

2. ¿Tendría sentido interpolar con parábolas?

Está claro que la interpolación lineal a trozos está abocada a ge-
nerar esquinas siempre que la tabla de datos no corresponda a una
función lineal. En general, la interpolación busca no solo una fun-
ción que pase por todos los puntos, sino que sea razonablemente suave
(por motivos no solo visuales sino de aproximación a la realidad).
Se podría intentar mejorar la aproximación lineal con funciones de
grado mayor, exigiendo que las tangentes en los puntos intermedios
fueran iguales. Podría intentarse con segmentos parabólicos: puesto
que tienen tres grados de libertad, al segmento i-ésimo se le puede
exigir que pase por los puntos ( xi−1 , yi−1 ) y ( xi , yi ) y que tenga la
misma derivada en xi que el siguiente. Esto puede parecer razonable
pero tiene cierta pega difícilmente salvable: genera una asimetría in-
trínseca al método (si uno plantea las ecuaciones de dicho sistema,
 3. SPLINES CÚBICOS: CURVATURA CONTINUA 69

f (x)
spline cuadrático
4

4 6 8 10

FIGURA 2. Interpolación cuadrática de una nube de 9

puntos. Compárese con la función original (en verde).

falta exactamente una por imponer para que sea compatible determi-
nado; esto hace que sea diferente el caso de un número par de nodos
y un número impar, o que la función de interpolación sea asimétri-
ca para datos simétricos). Sin entrar en detalles, puede comprobarse
que este spline cuadrático no es optimo, aunque aproxime la nube de
puntos de manera que no haya esquinas.
Tiene más problemas (por ejemplo, la curva se desvía mucho de la
nube de puntos si hay puntos cercanos en x pero lejanos en y en sen-
tidos diferentes). No se utiliza prácticamente nunca, salvo en cons-
trucción, por ejemplo, pues los arcos muy tendidos se aproximan con
parábolas.
El caso siguiente, el spline cúbico es el más utilizado: de hecho, es
lo que los programas de dibujo vectorial utilizan para trazar arcos
(aunque no con una tabla como (4), sino con dos tablas, pues las
curvas son parametrizadas como ( x (t), y(t))).

3. Splines cúbicos: curvatura continua

Para aproximar la nube de puntos con polinomios de grado tres,
se utiliza la siguiente definición:
DEFINICIÓN 14. Un spline cúbico de la tabla (4) es una función f :
[ x0 , xn ] → R tal que:
• La función f es un polinomio de grado 3 en cada segmento
[ xi−1 , xi ], para i = 1, . . . , n.
70 4. INTERPOLACIÓN

f (x)
spline cúbico
4

4 6 8 10

FIGURA 3. Spline cúbico de una nube de 9 puntos. Com-

párese con la función original (en verde).

• La función f es derivable dos veces con continuidad en todos

los puntos xi , para i = 1, . . . , n − 1.
De aquí se deduce que, si se llama Pi al polinomio de grado 3 que
coincide con f en [ xi−1 , xi ], entonces Pi′ ( xi ) = Pi′+1 ( xi ) y Pi′′ ( xi ) =
Pi′′+1 ( xi ); esto junto con el hecho de que Pi ( xi ) = yi y Pi ( xi+1 ) = yi+1
impone 4 condiciones para cada polinomio Pi , salvo para el primero
y el último, para los que solo hay 3 (esta es la simetría que poseen los
splines cúbicos y no poseen los cuadráticos). Por tanto, las condicio-
nes de spline cúbico determinan casi unívocamente los polinomios
Pi . Hace falta una condición en cada polinomio extremal (en P1 y
Pn ) para determinarlo totalmente. Estas condiciones pueden ser, por
ejemplo:
• Que la derivada segunda en los puntos extremos sea 0. A
esto se le denomina el spline cúbico natural, pero no tiene por
qué ser el mejor para un problema concreto. Las ecuaciones
son P1′′ ( x0 ) = 0 y Pn′′ ( xn ) = 0.
• Que la derivada tercera coincida en los puntos casi-extremos:
P1′′′ ( x1 ) = P2′′′ ( x1 ) y Pn′′′ ( xn−1 ) = Pn′′′−1 ( xn−1 ). Se dice que
este spline es extrapolado.
• Que haya cierta condición de periodicidad: P1′ ( x0 ) = Pn′ ( xn )
y lo mismo con la derivada segunda: P1′′ ( x0 ) = Pn′′ ( xn ). Es-
to tiene su interés, por ejemplo, si se está interpolando una
función periódica.
 3. SPLINES CÚBICOS: CURVATURA CONTINUA 71

3.1. El cálculo del spline cúbico: una matriz tridiagonal. Para

calcular efectivamente el spline cúbico conviene normalizar la nota-
ción. Como arriba, se llamará Pi al polinomio correspondiente al seg-
mento [ xi−1 , xi ] y se escribirá relativo al punto xi−1 :
Pi ( x ) = ai + bi ( x − xi−1 ) + ci ( x − xi−1 )2 + di ( x − xi−1 )3 ,
de modo que se han de calcular los ai , bi , ci y di a partir de la tabla de
datos (x, y) que se supone dada como arriba (4). La siguiente nor-
malización consiste en, en lugar de utilizar continuamente xi − xi−1 ,
llamar
hi = xi − xi−1 , para i = 1, . . . , n
(es decir, utilizar la anchura de los n intervalos en lugar de las coor-
denadas xi propiamente dichas).
Para aclarar todo el razonamiento, los datos conocidos se escribi-
rán en negrita.
El razonamiento es como sigue:
• Los ai se calculan directamente, pues Pi (xi−1 ) = ai y tiene
que ser igual a yi−1 . Por tanto,
ai = yi−1 para i = 1, . . . , n
• Se ha de cumplir que Pi (xi ) = yi , así que, utilizando la con-
dición anterior y llevando yi−1 al miembro de la derecha,
queda
(5) bi h i + c i h i 2 + d i h i 3 = y i − y i − 1 .
para i = 1, . . . , n. Esto da n ecuaciones.
• La condición de derivada continua es Pi′ (xi ) = Pi′+1 (xi ), así
que
(6) bi + 2ci hi + 3di hi 2 = bi+1 ,
para i = 1, . . . n − 1. De aquí salen n − 1 ecuaciones.
• Finalmente, las derivadas segundas han de coincidir en los
puntos intermedios, así que
(7) 2ci + 6di hi = ci+1 ,
para i = 1, . . . , n − 1. Esto da n − 1 ecuaciones.
En total se obtienen (aparte de las ai , que son directas), 3n − 2 ecua-
ciones para 3n incógnitas (las, b, c y d). Como ya se dijo, se imponen
condiciones en los extremos, pero al finalizar todo el razonamiento.
Una vez enunciadas las ecuaciones, se realizan simplificaciones
y sustituciones para obtener un sistema más inteligible. De hecho, se
72 4. INTERPOLACIÓN

despejan las d y las b en función de las c y queda un sistema lineal en

las c. Se procede como sigue:
Primero, se utiliza la ecuación (7) para despejar las di :
c − ci
(8) d i = i +1
3hi
y sustituyendo en (5), queda (hasta despejar bi ):
yi − yi−1 c + 2ci
(9) bi = − h i i +1 ;
hi 3
por otro lado, sustituyendo el di en (5) y operando hasta despejar bi ,
queda
(10) bi = bi − 1 + h i − 1 ( c i + c i − 1 ) .
para i = 2, . . . , n. Ahora no hay más que utilizar la ecuación (9) para
i y para i − 1 e introducirla en esta última, de manera que solo queden
c. Tras una colección de operaciones simples, se obtiene que
(11) ( )
yi − yi−1 yi−1 − yi−2
hi−1 ci−1 + (2hi−1 + 2hi )ci + hi ci+1 = 3 −
hi hi−1
para i = 2, . . . , n.
En forma matricial (ya sin utilizar negritas, para no recargar),
queda un sistema de ecuaciones Ac = α, donde A es
(12)
 
h1 2( h1 + h2 ) h2 0 ... 0 0 0
0 h2 2( h2 + h3 ) h3 . . . 0 0 0
A= .
 .. . . . . . .. 
.
.. .. .. . .. .. ..
0 0 0 0 . . . h n −1 2 ( h n −1 + h n ) h n
y c indica el vector columna de c1 , . . . , cn , mientras que α es el vector
columna  
α2
 α3 
 . 
 .. 
αn
con ( )
y i − y i −1 y i −1 − y i −2
αi = 3 − .
hi h i −1
Como se ve, este sistema es compatible indeterminado (le faltan dos
ecuaciones para tener solución única). Las ecuaciones se suelen im-
poner, como se explicó, como condiciones en los extremos. Por ejem-
plo, el spline natural significa que c1 = 0 y que cn = 0, así que A se
 3. SPLINES CÚBICOS: CURVATURA CONTINUA 73

completaría por arriba con una fila (1 0 0 . . . 0) y por abajo con otra
(0 0 0 . . . 1), mientras que al vector α se le añadiría una primera
componente igual a 0 y una última también igual a 0. Con estas n
ecuaciones se calcularían las ci y a partir de ellas, utilizando las ecua-
ciones (8) y (9), se calculan los valores de las demás variables.
El sistema de ecuaciones de un spline cúbico, somo se ve por la
ecuación (12), es tridiagonal: esto significa que solo tiene elementos
diferentes de cero en la diagonal principal y en las dos diagonales
adyacentes. (Esto es cierto para el spline natural y para cualquiera
que imponga condiciones en c1 y cn directamente). Estos sistemas se
resuelven muy fácilmente con una factorización LU —o bien, puede
implementarse directamente la solución como un algoritmo en fun-
ción de las αi y las hi . En cualquier caso, este tipo de sistemas tridia-
gonales es importantes reconocerlos y saber que su factorización LU
es rapidísima (pues si se piensa en el algoritmo de Gauss, cada ope-
ración de hacer ceros solo requiere hacerlo en una fila por debajo de la
diagonal).

3.2. El algoritmo: sencillo, resolviendo sistemas lineales. Tras

todo lo dicho, puede enunciarse el algoritmo para calcular el spline
cúbico que interpola una tabla de datos x, y de longitud n + 1, x =
( x0 , . . . , xn ), y = (y0 , . . . , yn ), en la que xi < xi+1 (y por tanto todos
los valores de x son diferentes), como indica el algoritmo 9.

3.3. Una cota sobre la aproximación. El hecho de que el spline

cúbico sea gráficamente satisfactorio no significa que sea técnicamente
útil. En realidad, lo es más de lo que parece. Si una función se “com-
porta bien” en la cuarta derivada, entonces el spline cúbico la apro-
xima razonablemente bien (y mejor cuanto más estrechos sean los
intervalos entre nodos). En concreto:

TEOREMA 6. Sea f : [ a, b] → R una función derivable 4 veces con

continuidad y supóngase que | f 4) ( x )| < M para x ∈ [ a, b]. Sea h el máximo
de las anchuras xi − xi−1 para i = 1, . . . , n. Si s( x ) es el spline cúbico para
los puntos ( xi , f ( xi )) que satisface que f ′ ( x0 ) = s′ ( x0 ) y f ′ ( xn ) = s′ ( xn ),
entonces
5M 4
|s( x ) − f ( x )| ≤ h .
384
Este resultado puede ser de gran utilidad, por ejemplo, para cal-
cular integrales o para acotar errores en soluciones de ecuaciones di-
ferenciales (que viene a ser lo mismo que integrar, vaya).
74 4. INTERPOLACIÓN

Algoritmo 9 (Cálculo del spline cúbico.)

Input: una tabla de datos x, y como se especifica arriba y dos con-
diciones sobre el primer nodo y el último (una en cada uno)
Output: un spline cúbico que interpola los puntos de dicha tabla.
Específicamente, una lista de n cuaternas ( ai , bi , ci , di ) tales que los
polinomios Pi ( x ) = ai + bi ( x − xi−1 ) + ci ( x − xi−1 )2 + di ( x − xi−1 )3
constituyen un spline cúbico para la tabla
⋆INICIO
ai ← yi−1 para i desde 1 hasta n
hi ← xi − xi−1 para i desde 1 hasta n
i←1
while i ≤ n do
if i > 1 and i < n then
Fi ← (0 · · · 0 hi−1 2(hi−1 + hi ) hi 0 · · · 0)
αi = 3(yi − yi−1 )/hi − 3(yi−1 − yi−2 )/hi−1
else
Fi ← la fila correspondiente a la ecuación para Pi
αi el coeficiente correspondiente a la condición para Pi
end if
i ← i+1
end while
M ← la matriz cuyas filas son las Fi para i = 1 hasta n
c ← M−1 α (resolver el sistema Mc = α)
i←1
while i < n do
bi ← (yi − yi−1 )/hi − hi (ci+1 + 2ci )/3
di ← (ci+1 − ci )/(3hi )
i ← i+1
end while
bn ← bn − 1 + h n − 1 ( c n + c n − 1 )
dn ← (yn − yn−1 − bn hn − cn h2n )/h3n
return ( a, b, c, d)

3.4. Definición de spline general. Se prometió incluir la defici-

nión general de spline:
DEFINICIÓN 15. Dada una nube de puntos como (4), un spline de
grado n, para n > 0, que los interpola es una función f : [ x0 , xn ] → R
tal que
• Pasa por todos los puntos: f ( xi ) = yi ,
• Es derivable n − 1 veces,
 4. INTERPOLACIÓN DE LAGRANGE 75

• En cada intervalo [ xi , xi+1 ] coincide con un polinomio de

grado n.
Es decir, es una función polinomial a trozos de grado n que pasa por
todos los puntos y es derivable n − 1 veces (se entiende que ser deri-
vable 0 veces significa ser continua).
De hecho, los que se utilizan son los de grado 1 y 3.

4. El polinomio interpolador de Lagrange: un solo polinomio

para todos los puntos
Hasta ahora se han estudiado splines, funciones que son polinó-
micas a trozos, para interpolar los datos de una tabla como (4), impo-
niendo la condición de que la función interpoladora pase por todos
los puntos de la tabla.
Se puede plantear el problema “límite”: buscar un polinomio que
pase por todos los puntos. Naturalmente, se procurará que sea de
grado mínimo (para, entre otras cosas, simplificar la búsqueda). Es
relativamente sencillo comprobar que existe uno de grado n (hay, re-
cuérdese, n + 1 datos) y que es único con esta condición. Este es el
polinomio interpolador de Lagrange:
TEOREMA 7 (del polinomio interpolador de Lagrange). Dada una
lista de puntos como (4) (recuérdese que xi < xi+1 ), existe un único poli-
nomio de grado n que pasa por cada ( xi , yi ) para i = 0, . . . , n.
La prueba del resultado es relativamente elemental. Tómese un
problema algo más “sencillo”: dados n + 1 valores x0 < · · · < xn , ¿se
puede construir un polinomio pi ( x ) de grado n que valga 1 en xi y ce-
ro en x j para j ̸= i? Con estas condiciones, se busca un polinomio de
grado n que tenga n ceros ya prefijados; es obvio que ha de ser múl-
tiplo de ϕi ( x ) = ( x − x0 )( x − x1 ) . . . ( x − xi−1 )( x − xi+1 ) . . . ( x − xn ).
Lo único que se ha de hacer es multiplicar a ϕi ( x ) por una constante
para que valga 1 en xi . El valor de ϕi ( x ) en xi es
ϕi ( xi ) = ( xi − x1 ) . . . ( xi − xi−1 )( xi − xi+1 ) . . . ( xi − xn ) = ∏ ( x i − x j ),
j ̸ =i

así que se ha de tomar

∏ j ̸ =i ( x − x j )
(13) pi ( x ) = .
∏ j ̸ =i ( xi − x j )
Estos polinomios, para i = 0, . . . , n, se llaman polinomios de base. Co-
mo se ha visto, cada uno de ellos vale 1 en el xi correspondiente y 0
76 4. INTERPOLACIÓN

en el resto, así que pueden pensarse como los vectores de la base es-
tándar de Rn+1 . El vector (y0 , y1 , . . . , yn ) es el que se quiere expresar
como combinación lineal de estos vectores, así que el polinomio que
pasa por los todos puntos ( xi , yi ) para i = 0, . . . , n es:
n
P ( x ) = y0 p0 ( x ) + y1 p1 ( x ) + · · · + y n p n ( x ) = ∑ yi pi ( x )
i =0
(14) n ∏ j ̸ =i ( x − x j )
= ∑ yi ∏ .
i =0 j ̸ =i ( xi − x j )

Para comprobar que es único, basta tener en cuenta que, si hubiera

otro, la diferencia sería un polinomio de grado n con n + 1 ceros, así
que la diferencia sería el polinomio nulo (y por tanto los dos que
pasan por ( xi , yi ) serían iguales).
El polinomio interpolador de Lagrange tiene algunas ventajas pe-
ro tiene un par de graves inconvenientes:
• Los denominadores que aparecen pueden ser muy peque-
ños cuando hay muchos puntos y dar lugar a errores de re-
dondeo.
• Es demasiado curvo.
El primer problema es intrínseco, pues los denominadores que
aparecen hay que calcularlos. El segundo depende esencialmente de
la distribución de puntos. Un ejemplo famoso e importante se debe a
Runge y ejemplifica el fenómeno de Runge: si se utilizan puntos equi-
distantes para interpolar una función con derivadas grandes en va-
lor absoluto, el polinomio interpolador de Lagrange se desvía mucho
(mucho, realmente) de la función, aunque pase por todos los puntos
de interpolación. Esto no le ocurre a los splines cúbicos (un spline
cúbico de once puntos de la función de Runge es indistinguible de
ella, por ejemplo).
Para solventar el problema de la curvatura del polinomio interpo-
lador de Lagrange cuando se trata de interpolar una función cono-
cida, se utilizan métodos basados en la idea de minimizar el valor
máximo que pueda tomar el polinomio

P( x ) = ( x − x0 )( x − x1 ) . . . ( x − xn ),

es decir, se busca una distribución de los puntos xi que resuelva un

problema de tipo minimax (minimizar un máximo). Sin entrar en de-
talles técnicos, dado el intervalo [−1, 1], se tiene que
 5. INTERPOLACIÓN APROXIMADA 77

0.8

0.6

0.4

0.2
1
f ( x ) = 1+12x 2
polin. interp.
0
−1 −0.5 0 0.5 1

FIGURA 4. Fenómeno de Runge: el polinomio interpo-

lador de Lagrange (rojo) se aleja ilimitadamente de la
función si los nodos están equiespaciados. Los trazos
negros siguen el spline cúbico (indistinguible de f ( x ).)

LEMA 9. Los puntos x0 , . . . , xn que minimizan el valor absoluto máximo

del polinomio P( x ) en el intervalo [−1, 1] vienen dados por la fórmula
( )
2k + 1 π
xi = cos
n+1 2
Por tanto, los puntos correspondientes para el intervalo [ a, b], donde a, b ∈
R, son
( b − a ) xi + ( a + b )
x̃i = .
2
A los puntos del lemma se les denomina nodos de Chebychev de un
intervalo [ a, b]. Son los que se han de utilizar si se quiere aproximar
una función por medio del polinomio interpolador de Lagrange.

5. Interpolación aproximada
En problemas relacionados con estudios estadísticos o experimen-
tales, los datos se suponen sujetos a errores, así que tratar de interpo-
lar una nube de puntos por medio de funciones que pasan por todos
ellos puede no tener ningún interes (de hecho, casi nunca lo tiene, en
este contexto). Para delimitar el problema con precisión, hace falta
indicar qué significa interpolar, pues cada problema puede tener una
idea distinta de lo que significa que una función se parezca a otra. En
el caso discreto (en el que se tiene una tabla de datos como (4), lo
habitual es intentar minimizar la distancia cuadrática de f ( xi ) a yi ,
78 4. INTERPOLACIÓN

0.8

0.6

0.4

0.2 1
f ( x ) = 1+12x 2
polin. interp.
0
−1 −0.5 0 0.5 1

FIGURA 5. Aproximación de la función de Runge por el

polinomio de Lagrange utilizando los nodos de Cheby-
chev. El error máximo cometido es mucho menor.

donde f sería la función de interpolación. Pero podría interesar algo

diferente (p.ej. minimizar la distancia de la gráfica de f a los puntos
( xi , yi ), un problema diferente del anterior), o cualquier otro criterio
—el que mejor corresponda a los datos.

5.1. Interpolación por mínimos cuadrados. El método más uti-

lizado de interpolación, sin duda alguna, es el de mínimos cuadrados.
Se parte de una nube de puntos x, y, donde x e y son vectores de ta-
maño n arbitrarios (es decir, puede haber valores repetidos en las x
y el orden es irrelevante). Dada una función f : R → R, se define:
DEFINICIÓN. El error cuadrático de f en el punto xi (una componen-
te de x) es el valor ( f ( x ) − yi )2 . El error cuadrático total de f en la nube
dada por x e y es la suma
n
∑ ( f ( x i ) − y i )2 .
i =0

El problema de la interpolación por mínimos cuadrados consiste

en, dada la nube de puntos, encontrar, en un conjunto determinado de
funciones una que minimice el error cuadrático total. El matiz de “en
un conjunto de funciones” es crucial. De hecho, en estas notas se estu-
diará con precisión el problema en un espacio vectorial de funciones
y se tratará de manera aproximada el de algunos conjuntos que no
son espacios vectoriales.
 5. INTERPOLACIÓN APROXIMADA 79

1
−1 −0.5 0 0.5 1

FIGURA 6. Interpolación de una nube de puntos por mí-

nimos cuadrados lineales por una ecuación del tipo
y = ax + b.

5.1.1. Mínimos cuadrados lineales. Supóngase que ha de aproximar-

se una nube de puntos de tamaño N mediante una función f que
pertenece a un espacio vectorial V de funciones y que se conoce una
base de V: f 1 , . . . , f n . Se supone siempre que N > n (y de hecho, ha-
bitualmente, N es mucho mayor que n). Es decir, se busca la función
f ∈ V tal que
N
∑ ( f ( x i ) − y i )2
i =1

es mínimo. Pero, puesto que V es un espacio vectorial, dicha función

es una combinación lineal de los elementos de la base:
f = a1 f 1 + a2 f 2 + · · · + a n f n .
Y, en realidad, se puede plantear el problema como la búsqueda de
los coeficientes a1 , . . . , an que hacen mínima la función
N
F ( a1 , . . . , a n ) = ∑ ( a1 f 1 ( x i ) + · · · + a n f n ( x i ) − y i )2 ,
i =1

es decir, dada F ( a1 , . . . , an ), una función de n variables, se ha de cal-

cular un mínimo. La expresión de F indica claramente que es una
función diferenciable; por tanto, el punto que dé el mínimo resuelve
las ecuaciones correspondientes a hacer las parciales iguales a 0. Es
80 4. INTERPOLACIÓN

decir, hay que resolver el sistema

∂F ∂F ∂F
= 0, = 0, . . . , = 0.
∂a1 ∂a2 ∂an
Si se escribe la parcial de F con respecto a a j , queda
N
∂F
= ∑ 2 f j ( xi )( a1 f 1 ( xi ) + · · · + an f n ( xi ) − yi ) = 0.
∂a j i =1

Si, para simplificar, se denomina

N
yj = ∑ f j ( xi ) yi ,
i =1

uniendo todas las ecuaciones y escribiendo en forma matricial, sale

el sistema
   
f 1 ( x1 ) f 1 ( x2 ) ... f1 (x N ) f 1 ( x1 ) f 2 ( x1 ) ... f n ( x1 ) a1
 f 2 ( x1 ) f 2 ( x2 ) ... f2 (x N )   f 1 ( x2 ) f 2 ( x2 ) ... f n ( x2 )   a2 
   
 . .. .. ..   .. .. .. ..   ..  =
 .. . . .  . . . .  . 
f n ( x1 ) f n ( x2 ) ... fn (xN ) f1 (x N ) f2 (x N ) ... fn (xN ) an
(15)  
y1
 y2 
 
=  . ,
 .. 
yn

que puede escribirse

XX t a = Xy
donde se sobreentiende que la matriz X es la lista (por filas) de los va-
lores de cada f i en los puntos x j , y la y es el vector columna (y0 , y1 , . . . , yn ) T ,
es decir:
   
f 1 ( x1 ) f 1 ( x2 ) . . . f 1 ( x N ) y0
 f 2 ( x1 ) f 2 ( x2 ) . . . f 2 ( x N )   y1 
X=  ... .. .. ..  
.
, y= .. 
. . .
f n ( x1 ) f n ( x2 ) . . . fn (xN ) yn

Este sistema es compatible determinado si hay al menos tantos pun-

tos como la dimensión de V y las funciones f i generan filas indepen-
dientes en la matriz X.
Es bastante probable que el sistema dado por (15) no esté muy
bien acondicionado.
 6. CÓDIGO DE ALGUNOS ALGORITMOS 81

5.1.2. Interpolación no lineal: aproximaciones. En bastantes casos se

requiere aproximar una nube de puntos por una función que perte-
necza a una familia que no forme un espacio vectorial. Un ejemplo
clásico es tratar de encontrar una función del tipo
2
(16) f ( x ) = aebx
que aproxime una cierta nube de puntos. Las funciones de ese tipo
dan lugar (para a, b > 0) a “campanas de Gauss”. Está claro que di-
chas funciones no forman un espacio vectorial, así que no se pueden
aplicar las técnicas del apartado anterior.
Sin embargo, puede intentar trasladarse el problema a una familia
lineal, aproximar esta por mínimos cuadrados y calcular la función
original equivalente.
Por seguir con el ejemplo, si se toman logaritmos en ambos lados
de la ecuación (16), se obtiene la expresión
log( f ( x )) = log( a) + bx2 = a′ + bx2 ,
y, considerando las funciones 1 y x2 , se puede ahora aproximar la
nube de puntos x, log(y), en lugar de la original, utilizando mínimos
cuadrados lineales (pues 1 y x2 forman un espacio vectorial). Si se
obtiene la solución log( a0 ), log(b0 ), entonces puede pensarse que
2
g( x ) = e a0 eb0 x
será una buena aproximación de la nube de puntos original x, y. Pe-
ro es importante ser consciente de que esta aproximación posiblemente
no sea la mejor respecto del error cuadrático total. Piénsese que, al
tomar logaritmos, los puntos de la nube cercanos a 0 estarán cerca
de −∞, con lo cual tendrán un peso mayor que los que al principio
están cerca de 1 (que pasan a estar cerca de 0): al tomar logaritmos,
los errores absolutos y relativos cambian de manera no lineal. Más
aun, si uno de los valores de x es 0, el problema no puede convertirse
con el logaritmo, y ha de eliminarse ese valor o cambiarse por otro
aproximado…

6. Código de algunos algoritmos

6.1. El spline cúbico natural. Matlab, por defecto, calcula spli-
nes cúbicos con la condición not-a-knot, que es cierta relación entre
las derivadas terceras en los puntos segundo y penúltimo. Solo se
pueden calcular splines naturales con un toolbox. El código que sigue
implementa el cálculo del spline cúbico natural de una colección de
puntos. La entrada es:
82 4. INTERPOLACIÓN

3 datos
3e−2x
2

0
−2 −1 0 1 2

FIGURA 7. Interpolación no lineal tomando logaritmos:

la nube de puntos se parece mucho a la función f ( x ),
pero la interpolación lineal tomando logaritmos y lue-
go haciendo la exponencial es muy mala (verde). El pro-
blema reside en los valores de y muy cercanos a 0: al
tomar logaritmos, adquieren una importancia despro-
porcionada.

x: la lista de las coordenadas x de los puntos

y: la lista de las coordenadas y de los puntos
Devuelve un “objeto” de Matlab que se denomina polinomio a tro-
zos, que describe exactamente un objeto polinomial definido a trozos:
los intervalos en los que está definido vienen dados por el vector x y
su valor en un t (que hay que computar utilizando la función ppval)
viene dado por el valor del polinomio correspondiente al intervalo
de las x que contiene a t.
Un ejemplo de uso podría ser el siguiente, para comparar el spline
cúbico con la gráfica de la función seno:
> x = linspace(-pi , pi , 10);
> y = sin(x);
> f = spline_cubico(x, y);
> u = linspace(-pi , pi , 400);
> plot(u, ppval(f, u)); hold on;
> plot(u, sin(u), 'r');

6.2. El polinomio interpolador de Lagrange. El cálculo del Po-

linomio interpolador de Lagrange en Matlab/Octave es extremada-
mente simple, como se puede ver mirando el Código 9.
La entrada es la siguiente:
 6. CÓDIGO DE ALGUNOS ALGORITMOS 83

% natural cubic spline: second derivative at both

% endpoints is 0. Input is a pair of lists describing
% the cloud of points.
function [f] = spline_cubico(x, y)
% safety checks
n = length(x) -1;
if(n<=1 | length(x) ~= length(y))
8 warning('Wrong data ')
f= [];
return
end

% variables and coefficients for the linear system ,

% these are the ordinary names. Initialization
a = y(1:n);
h = diff(x);
F = zeros(n);
18 alpha = zeros(n, 1);

% true coefficients (and independent terms) of the linear system

for k=1:n
if(k> 1& k < n)
F(k,[k-1 k k+1]) = [h(k-1), 2*(h(k-1) + h(k)), h(k)] ;
alpha(k) = 3*(y(k+1)-y(k))/h(k) - 3*(y(k) - y(k-1))/h(k-1);
else
% these two special cases are the 'natural ' condition
% (second derivatives at endpoints = 0)
28 F(k,k) = 1;
alpha(k) = 0;
end
k=k+1;
end

% These are the other coefficients of the polynomials

% (a + bx + cx^2 + dx^3) ... Initialization
c = (F\alpha) ';
b = zeros(1,n);
38 d = zeros(1,n);
% unroll all the coefficients as in the theory
k = 1;
while(k<n)
b(k) = (y(k+1)-y(k))/h(k) - h(k) *(c(k+1) +2*c(k))/3;
k=k+1;
end
d(1:n-1) = diff(c)./(3*h(1:n-1));

% the last b and d have explicit expressions:

48 b(n) = b(n-1) + h(n-1)*(c(n)+c(n-1));
d(n) = (y(n+1)-y(n)-b(n)*h(n)-c(n)*h(n)^2)/h(n)^3;

% finally , build the piecewise polynomial (a Matlab function)

% we might implement it by hand , though
f = mkpp(x,[d; c; b ;a ]');
end

FIGURA 8. Código del algoritmo para calcular un spline cúbico

 84 4. INTERPOLACIÓN

% Lagrange interpolation polynomial

% A single base polynomial is computed at
% each step and then added (multiplied by
4 % its coefficient) to the final result.
% input is a vector of x coordinates and
% a vector (of equal length) of y coordinates
% output is a polynomial in vector form (help poly).
function [l] = lagrange(x,y)
n = length(x);
l = 0;
for m=1:n
b = poly(x([1:m-1 m+1:n]));
c = prod(x(m)-x([1:m-1 m+1:n]));
14 l = l + y(m) * b/c;
end
end

FIGURA 9. Código del algoritmo para calcular el polino-

mio interpolador de Lagrange

x: El vector de las coordenadas x de la nube de puntos que se

quiere interpolar
y: El vector de las coordenadas y de la nube
La salida es un polinomio en forma vectorial, es decir, una lista de
los coeficientes an , an−1 , . . . , a0 tal que el polinomio P( x ) = an x n +
· · · + a1 x + a0 es el polinomio interpolador de Lagrange para la nube
de puntos (x, y).

6.3. Interpolación lineal. Para la interpolación lineal en Matla-

b/Octave se ha de utilizar un objeto especial: un array de celdas: la
lista de funciones que forman la base del espacio vectorial se ha de
introducir entre llaves.
La entrada es la siguiente:
x: El vector de las coordenadas x de la nube de puntos que se
quiere interpolar
y: El vector de las coordenadas y de la nube
F: Un array de celdas de funciones anónimas. Cada elemento
es una función de la base del espacio vectorial que se utiliza
para interpolar
La salida es un vector de coeficientes c, tal que c1 F1 + · · · + cn Fn es
la función de interpolación lineal de mínimos cuadrados de la nube
(x, y).
Un ejemplo de uso podría ser el siguiente, para aproximar con
funciones trigonométricas
 6. CÓDIGO DE ALGUNOS ALGORITMOS 85

% interpol.m
% Linear interpolation.
% Given a cloud (x,y), and a Cell Array of functions F,
4 % return the coefficients of the least squares linear
% interpolation of (x,y) with the base F.
#
% Input:
% x: vector of scalars
% y: vector of scalars
% F: Cell array of anonymous functions
#
% Outuput:
% c: coefficients such that
14 % c(1)*F{1,1} + c(2)*F{1,2} + ... + c(n)*F{1,n}
% is the LSI function in the linear space <F>.

function [c] = interpol(x, y, F)

n = length(F);
m = length(x);
X = zeros(n, m);
for k=1:n
X(k,:) = F{1,k}(x);
end
24 A = X*X.';
b = X*y.';
c = (A\b) ';
end

FIGURA 10. Código del algoritmo para calcular la inter-

polación lineal por mínimos cuadrados.

> f1=@(x) sin(x); f2=@(x) cos(x); f3=@(x) sin (2*x); f4=@(x) cos (2*x);
> f5=@(x)sin (3*x); f6=@(x)cos (3*x);
> F={f1 , f2 , f3 , f4 , f5 , f6};
> u=[1,2,3,4,5,6];
> r=@(x) 2*sin(x)+3* cos(x) -4*sin (2*x) -5*cos (3*x);
> v=r(u)+rand (1,6) *.01;
> interpol(u,v,F)
ans =
1.998522 2.987153 -4.013306 -0.014984 -0.052338 -5.030067
 CAPÍTULO 5

Derivación e Integración Numéricas

Se expone con brevedad la aproximación de la derivada de una

función mediante la fórmula simétrica de incrementos y cómo este
tipo de fórmula puede generalizarse para derivadas de órdenes su-
periores. El problema de la inestabilidad de la derivación numérica
se explica someramente.
La integración se trata con más detalle (pues el problema es mu-
cho más estable que el anterior), aunque también con brevedad. Se
definen las nociones de fórmulas de cuadratura de tipo interpolato-
rio y de grado de precisión, y se explican las más sencillas (trapecios
y Simpson) en el caso simple y en el compuesto.

1. Derivación Numérica: un problema inestable

A veces se requiere calcular la derivada aproximada de una fun-
ción mediante fórmulas que no impliquen la derivada, bien sea por-
que esta es difícil de computar —costosa—, bien porque realmente
no se conoce la expresión simbólica de la función a derivar. En el pri-
mer caso se hace necesario sobre todo utilizar fórmulas de derivación
numérica lo más precisas posibles, mientras que en el segundo caso
lo que se requerirá será, casi con certeza, conocer una cota del error
cometido. En estas notas solo se va a mostrar cómo la fórmula simétri-
ca para calcular derivadas aproximadas es mejor que la aproximación
naíf consistente en calcular el cociente del incremento asimétrico. Se
expone la fórmula equivalente para la derivada segunda y se da una
idea de la inestabilidad del método.

1.1. La derivada aproximada simétrica. Una idea simple para

calcular la derivada de una función f en un punto x es, utilizando
la noción de derivada como límite de las secantes, tomar la siguiente
aproximación:
f ( x + h) − f ( x )
(17) f ′ (x) ≃ ,
h
donde h será un incremento pequeño de la variable independiente.
87
88 5. DERIVACIÓN E INTEGRACIÓN NUMÉRICAS

f ( x0 + h )

f ( x0 )

f ( x0 − h )

x0 − h x0 x0 + h

FIGURA 1. Derivada aproximada: por la derecha (rojo),

por la izquierda (verde) y centrada (amarillo). La cen-
trada es mucho más parecida a la tangente (magenta).

Sin embargo, la propia fórmula (17) muestra su debilidad: ¿se ha

de tomar h positivo o negativo? Esto no es irrelevante. Supongamos
que f ( x ) = 1/x y se quiere calcular la derivada en x = 1/2. Utilice-
mos aproximaciones de |h| = .01. Sabemos que f ′ (0.5) = 0.25. Con
la aproximación “natural”, se tiene, por un lado,
f ( x + .01) − f ( x ) 1
− 12
= 2.01
= −0,248756+
.01 .01
mientras que, por otro,
f ( x − .01) − f ( x ) 1
− 1
= − 1.99 2
= −0,251256+
−.01 .01
que difieren en el segundo decimal. ¿Hay alguna de las dos que tenga
preferencia en un sentido abstracto? No.
Cuando se tienen dos datos que aproximan un tercero y no hay
ninguna razón abstracta para elegir uno por encima de otro, parece
razonable utilizar la media aritmética. Probemos en este caso:
( )
1 f ( x + h) − f ( x ) f ( x − h) − f ( x )
+ = −0.2500062+
2 h −h
que aproxima la derivada real 0.25 con cuatro cifras significativas
(cinco si se trunca).
Esta es, de hecho, la manera correcta de plantear el problema de
la derivación numérica: utilizar la diferencia simétrica alrededor del
 2. INTEGRACIÓN NUMÉRICA: FÓRMULAS DE CUADRATURA 89

punto y dividir por dos veces el intervalo (es decir, la media entre la
“derivada por la derecha” y la “derivada por la izquierda”).
TEOREMA 8. La aproximación naíf a la derivada tiene una precisión de
orden 1, mientras que la fórmula simétrica tiene una precisión de orden 2.
PRUEBA. Sea f una función con derivada tercera en el intervalo
[ x − h, x + h]. Utilizando el Polinomio de Taylor de grado uno, se tiene
que, para cierto ξ ∈ [ x − h, x + h],
f ′′ (ξ ) 2
f ( x + h) = f ( x ) + f ′ ( x )h + h ,
2
así que, despejando f ′ ( x ), queda
f ( x + h) − f ( x ) f ′′ (ξ )
f ′ (x) = − h,
h 2
que es lo que significa que la aproximación tiene orden 2. Como se
ve, el término de más a la derecha no puede desaparecer.
Sin embargo, para la fórmula simétrica se utiliza el Polinomio de
Taylor dos veces, y de segundo grado:
f ′′ ( x ) 2 f 3) ( ξ ) 3
f ( x + h) = f ( x ) + f ′ ( x )h + h + h ,
2 6
f ′′ ( x ) 2 f 3) ( ζ ) 3
f ( x − h) = f ( x ) − f ′ ( x )h + h − h
2 6
para ξ ∈ [ x − h, x + h] y ζ ∈ [ x − h, x + h]. Restando ambas igualda-
des queda
f ( x + h) − f ( x − h) = 2 f ′ ( x )h + K (ξ, ζ )h3
donde K es una mera suma de los coeficientes de grado 3 de la ecua-
ción anterior, así que
f ( x + h) − f ( x − h)
f ′ (x) = + K (ξ, ζ )h2 ,
2h
que es justo lo que significa que la fórmula simétrica es de segundo
orden. □

2. Integración Numérica: fórmulas de cuadratura

La integración numérica, como corresponde a un problema en el
que los errores se acumulan (un problema global, por así decir) es más
estable que la derivación, por extraño que parezca (pese a que integrar
simbólicamente es más difícil que derivar, la aproximación numérica
se comporta mucho mejor).
90 5. DERIVACIÓN E INTEGRACIÓN NUMÉRICAS

En general, uno está interesado en dar una fórmula de integración

numérica. Es decir, se busca una expresión algebraica general tal que,
dada una función f : [ a, b] → R cuya integral se quiere aproximar,
se pueda sustituir f en la expresión y obtener un valor —la integral
aproximada.
DEFINICIÓN 16. Una fórmula de cuadratura simple para un intervalo
[ a, b] es una colección de puntos x1 < x2 < · · · < xn+1 ∈ [ a, b] y
de valores a1 , . . . , an+1 . La integral aproximada de una función f en [ a, b]
utilizando dicha fórmula de cuadratura es la expresión
a 1 f ( x 1 ) + a 2 f ( x 2 ) + · · · + a n +1 f ( x n +1 ).
Es decir, una fórmula de cuadratura no es más que “la aproxima-
ción de una integral utilizando valores intermedios y pesos”. Se dice
que la fórmula de cuadratura es cerrada si x1 = a y xn+1 = b; si ni a
ni b están entre los xi , se dice que la fórmula es abierta.
Si los puntos xi están equiespaciados, la fórmula se dice que es de
Newton-Coates.
Lógicamente, interesa encontrar las fórmulas de cuadratura que
mejor aproximen integrales de funciones conocidas.
DEFINICIÓN 17. Se dice que una fórmula de cuadratura es de orden
n si para cualquier polinomio P( x ) de grado n se tiene que
∫ b
P( x ) dx = a1 P( x1 ) + a2 P( x2 ) + · · · + an+1 P( xn+1 ).
a
Es decir, si integra con exactitud polinomios de grado n.
Las fórmulas de cuadratura más sencillas son: la del punto medio
(abierta, n = 1), la fórmula del trapecio (cerrada, n = 2) y la fórmula
de Simpson (cerrada, n = 3). Se explican a continuación.
Se explican la fórmulas simples a continuación y luego se genera-
lizan a sus versiones compuestas.
2.1. La fórmula del punto medio. Una manera burda pero natu-
ral de aproximar una integral es multiplicar el valor de la función en
el punto medio por la anchura del intervalo. Esta es la fórmula del
punto medio:
DEFINICIÓN 18. La fórmula de cuadratura del punto medio corres-
ponde a x1 = ( a + b)/2 y a1 = (b − a). Es decir, se aproxima
∫ b ( )
a+b
f ( x ) dx ≃ (b − a) f .
a 2
por el área del rectángulo que tiene un lado horizontal a la altura del
valor de f en el punto medio.
 2. INTEGRACIÓN NUMÉRICA: FÓRMULAS DE CUADRATURA 91

f (x)
2

0
2 3 3.5 4 5

FIGURA 2. Fórmula del punto medio: se aproxima el

área debajo de f por el rectángulo.

Es elemental comprobar que la fórmula del punto medio es de orden

1: integra con precisión funciones lineales pero no polinomios de se-
gundo grado.
2.2. La fórmula del trapecio. La siguiente aproximación natural
(no necesariamente óptima) es utilizar dos puntos. Puesto que ya hay
dos dados, a y b, lo intuitivo es utilizarlos.
Dados dos puntos a y b, se tendrán los valores de f en cada uno
de ellos. Se puede interpolar linealmente f con una recta y aproximar
la integral por medio de esta o bien aproxmar la integral con dos rec-
tángulos como en la regla del punto medio y hacer la media. Ambas
operaciones producen el mismo resultado. En concreto:
DEFINICIÓN 19. La fórmula del trapecio para [ a, b] consiste en tomar
x1 = a, x2 = b y pesos a1 = a2 = (b − a)/2. Es decir, se aproxima
∫ b
b−a
f ( x ) dx ≃ ( f ( a) + f (b))
a 2
la integral de f por el área del trapecio con base en OX, que une
( a, f ( a)) con (b, f (b)) y es paralelo al eje OY.
La fórmula del trapecio, pese a incorporar un punto más que la
regla del punto medio, es también de orden 1.
2.3. La fórmula de Simpson. El siguiente paso natural (que no
óptimo) es utilizar, en lugar de 2 puntos, 3 y en vez de aproximar
mediante rectas, utilizar parábolas. Esta aproximación es singular-
mente eficaz, pues tiene orden 3. Se denomina fórmula de Simpson.
92 5. DERIVACIÓN E INTEGRACIÓN NUMÉRICAS

f (x)
2

0
2 3 4 5

FIGURA 3. Fórmula del trapecio: se aproxima el área de-

bajo de f por el trapecio.

DEFINICIÓN 20. La fórmula de Simpson es la fórmula de cuadratura

que consiste en tomar x1 = a, x2 = ( a + b)/2 y x3 = b, y como pesos,
los correspondientes a la interpolación correcta de un polinomio de
4( b − a )
grado 2. Es decir1, a1 = b−6 a , a2 = 6 y a3 = b−6 a . Así pues, consiste
en aproximar
∫ b ( ( ) )
b−a a+b
f ( x ) dx ≃ f ( a) + 4 f + f (b)
a 6 2
la integral por una media ponderada de áreas de rectángulos inter-
medios. Esta fórmula es sorprendentemente precisa y conviene sabér-
sela de memoria: un sexto de la anchura por la ponderación 1, 4, 1 de
los valores en extremo, punto medio, extremo.
Lo singular de la fórmula de Simpson es que es de tercer orden: pese
a interpolar con parábolas, se obtiene una fórmula que integra con pre-
cisión polinomios de tercer grado. Para terminar, téngase en cuenta
que la ecuación de la parábola que pasa por los tres puntos es irrelevante: no
hace falta calcular el polinomio interpolador, basta usar la fórmula de
la Definición 20.

2.4. Las fórmulas compuestas. Las fórmulas de cuadratura com-

puestas consisten en aplicar una fórmula de cuadratura en subinter-
valos. En lugar de aproximar, por ejemplo, utilizando la fórmula del

1Esto es sencillo de calcular, usando un polinomio de grado 3, por ejemplo.

 2. INTEGRACIÓN NUMÉRICA: FÓRMULAS DE CUADRATURA 93

f (x)
2

0
2 3 3.5 4 5

FIGURA 4. Fórmula de Simpson: se aproxima el área de-

bajo de f por la parábola que pasa por los puntos ex-
tremos y el punto medio (en negro).

trapecio
∫ b
b−a
f ( x ) dx ≃ ( f ( a) + f (b))
a 2
se subdivide el intervalo [ a, b] en subintervalos y se aplica la fórmula
en cada uno de estos.
Podría hacerse sin más, mécanicamente, pero el caso es que, si la
fórmula es cerrada, para las fórmulas de Newton-Coates, siempre es
más sencillo aplicar la fórmula general que ir subintervalo a subintervalo,
pues los valores en los extremos se ponderan.
Nota: aún así, al alumno se le aconseja no aprenderse las fórmulas
siguientes de memoria, pues lo que suele ocurrir es que se haga un
lío. Basta con aplicarlas subintervalo a subintervalo. Los siguientes
apartados pueden omitirse, y de hecho es lo que se aconseja. De ahí
que estén en gris claro.
2.4.1. Fórmula compuesta de los trapecios. Si se tienen dos intervalos
consecutivos [ a, b] y [b, c] de igual longitud (es decir, b − a = c − b) y se
aplica la fórmula de los trapecios en [ a, b] y en [b, c] para aproximar
la integral total de f entre a y c, queda:
∫ c
b−a c−b
f ( x ) dx = ( f (b) + f ( a)) + ( f (c) + f (b)) =
a 2 2
h
( f ( a) + 2 f (b) + f (c)) ,
2
94 5. DERIVACIÓN E INTEGRACIÓN NUMÉRICAS

f (x)
2

0
2 2.5 3 3.5 4 4.5 5

FIGURA 5. Fórmula del trapecio compuesta: simple-

mente repetir la fórmula del trapecio en cada subinter-
valo.

donde h = b − a es la anchura del subintervalo, que es una ponderación

de los valores en puntos extremos e intermedio). En el caso general,
queda algo totalmente análogo:
DEFINICIÓN 21 (Fórmula de los trapecios compuesta). Dado un
intervalo [ a, b] y un número de nodos n + 1, la fórmula de los trape-
cios compuesta para [ a, b] con n nodos viene dada por los nodos x0 =
a, x1 = a + h, . . . , xn = b, con h = (b − a)/(n − 1)2 y por la aproxi-
mación
∫ b
h
f ( x ) dx ≃ ( f ( a) + 2 f ( x1 ) + 2 f ( x2 ) + · · · + 2 f ( xn−1 ) + f (b)) .
a 2
Es decir, la mitad de la anchura de cada intervalo por la suma de
los valores en los extremos y los dobles de los valores en los puntos
interiores.
2.4.2. Fórmula de Simpson compuesta. De modo análogo, la fórmu-
la de Simpson compuesta consiste en utilizar la fórmula de Simpson
en varios intervalos que se unen para formar uno mayor. Como en la
anterior, al ser iguales los extremos final e inicial de intervalos con-
secutivos, la fórmula final compuesta no es una simple suma tonta de
todas las fórmulas intermedias: hay cierta combinación que simplifi-
ca el resultado final.
2Esto no es más que la anchura de cada subintervalo, no hay que liarse con la
n.
 2. INTEGRACIÓN NUMÉRICA: FÓRMULAS DE CUADRATURA 95

80 f (x)

60

40

20

0
2 2.75 3.5 4.25 5

FIGURA 6. Fórmula de Simpson compuesta: concep-

tualmente, repetir la fórmula de Simpson en cada
subintervalo.

DEFINICIÓN 22 (Fórmula de Simpson compuesta). Dado un inter-

valo [ a, b] dividido en n subintervalos (y por tanto, dado un número
de nodos equiespaciados 2n + 1), la fórmula de Simpson compuesta pa-
ra [ a, b] con 2n + 1 nodos viene dada por los nodos x0 = a, x1 =
a + h, . . . , x2n = b, (así que h = (b − a)/(2n)) y por la aproximación
∫ b
b−a
f ( x ) dx ≃ ( f ( a ) + 4 f ( x1 ) + 2 f ( x2 ) + 4 f ( x3 ) + . . .
a 6n
· · · + 2 f ( x2n−2 ) + 4 f ( x2n−1 ) + f (b)).
Es decir, un tercio del espacio entre nodos multiplicado por: el valor
en los extremos más el doble de la suma de los valores en los puntos
extremos de los subintervalos más cuatro veces la suma de los valores
en los puntos medios de los subintervalos.
Es decir, la fórmula de Simpson compuesta es la misma que la
de Simpson en todos los aspectos salvo en los puntos extremos de los
intervalos, donde hay que multiplicar por dos —obviamente, pues se
están sumando. Es más fácil recordar la fórmula si simplemente se
considera que hay que aplicar la fórmula simple en cada subintervalo
y sumar todos los términos que salen.
 CAPÍTULO 6

Ecuaciones diferenciales

Sin duda, la actividad primordial para la que se utilizan los méto-

dos numéricos es integrar ecuaciones diferenciales (es la denominación
clásica del proceso de resolución de una ecuación diferencial).

1. Introducción
Una ecuación diferencial es un tipo especial de ecuación funcio-
nal: una ecuación (es decir, una igualdad en la que alguno de los
elementos es desconocido) en la que una de las incógnitas es una
función. Ya se han resuelto ecuaciones diferenciales antes: cualquier
integral es la solución de una ecuación diferencial. Por ejemplo,
y′ = x
es una ecuación en la que se busca conocer una función y( x ) cuya
derivada es x. Como se sabe, la solución no es única: hay una constante
de integración, de manera que la solución general se escribe
x2
y= + C.
2
En realidad, la constante de integración tiene una manera más
natural de entenderse. Cuando se calcula una primitiva, lo que se está
buscando es una función cuya derivada es conocida. Una función, al
fin y al cabo, no es más que una gráfica en el plano X, Y. Esta gráfica
puede estar situada más arriba o más abajo (la derivada no cambia
por dibujar una función paralelamente arriba o abajo). Sin embargo,
si lo que se busca es resolver el problema
y′ = x, y(3) = 28,
entonces se está buscando una función cuya derivada es x con una
condición puntual: que su valor en 3 sea 28. Una vez que este valor
está fijado, solo hay una gráfica que tenga esa forma y pase por ese punto (el
(3, 28)): a esto se le llama una condición inicial: se busca una función
cuya derivada cumpla algo, pero imponiendo que pase por determinado
97
98 6. ECUACIONES DIFERENCIALES

punto: de este modo la constante ya no aparece, se puede determinar

sin más sustituyendo:
28 = 32 + C ⇒ C = 19.
Esta es la idea de condición inicial para una ecuación diferencial.
Considérese la ecuación
y′ = y
(cuya solución general debería ser obvia). Esta ecuación significa: la
función y( x ) cuya derivada es igual a ella misma. Uno tiende siempre
a pensar en y( x ) = e x , pero ¿es esta la única solución? Si el problema
se piensa de otra manera, geométrica, la ecuación significa: la fun-
ción y( x ) cuya derivada es igual a la altura (y) en cada punto de la
gráfica. Pensada así, está claro que tiene que haber más de una fun-
ción solución. De hecho, uno supone que por cada punto del plano
( x, y) pasará al menos una solución (pues geométricamente no hay
ninguna razón por la que los puntos ( x, e x ) sean preferibles a otros
para que pase por ellos una solución.) En efecto, la solución general
es
y( x ) = Ce x ,
donde C es una constante de integración. Ahora, si especificamos una
condición inicial, que y( x0 ) = y0 para x0 , y0 ∈ R, entonces está claro
que
y0 = Ce x0
y por tanto C = y0 /e x0 , que está siempre definido (pues el denomi-
nador no es 0.)
De este modo, dada una ecuación diferencial, hace falta especifi-
car al menos una condición inicial para que la solución esté completa-
mente determinada. Si se considera
y′′ = −y,
no basta con especificar una sola condición inicial. ¿Qué soluciones
tiene esta ecuación? En breve, son las de la forma y( x ) = a( x ) +
b cos( x ), para dos constantes a, b ∈ R. Para que la solución esté total-
mente determinada hay que fijar esas dos constantes. Esto, por ejem-
plo, puede hacerse imponiendo (es lo común) el valor en un punto,
y(0) = 1 y el valor de la derivada en ese punto, y′ (0) = −1. En cuyo
caso, la solución es y( x ) = −( x ) + cos( x ).
Así que este tema trata de la solución aproximada de ecuaciones
diferenciales. De momento se estudiarán solo funciones de una va-
riable y ecuaciones en las que solo aparece la primera derivada, pero
esto puede variar en el futuro.
 2. LO MÁS BÁSICO 99

2. Lo más básico
Comencemos por el principio:
DEFINICIÓN 23. Una ecuación diferencial ordinaria es una igualdad
a = b en la que aparece la derivada (o una derivada de orden supe-
rior, pero no derivadas parciales) de una variable libre cuyo rango es
el conjunto de funciones de una variable real.
Obviamente, así no se entiende nada. Se puede decir que una
ecuación diferencial ordinaria es una ecuación en la que una de las in-
cógnitas es una función de una variable y( x ), cuya derivada aparece
explícitamente en dicha ecuación. Las no ordinarias son las ecuaciones
en derivadas parciales, que no se tratan en este texto, de momento.
EJEMPLO 12. Arriba se han dado dos ejemplos, las ecuaciones di-
ferenciales pueden ser de muchos tipos:
y′ = sin( x )
xy = y′ − 1
(y′ )2 − 2y′′ + x2 y = 0
y′
− xy = 0
y
etc.
De ahora en adelante se supondrá que la ecuación diferencial de
partida solo tiene una variable libre, que es la de la función que se
busca, que se asume que es y.
DEFINICIÓN 24. Se dice que una ecuación diferencial tiene orden n
si n es la mayor derivada de y que aparece. Si y y sus derivadas solo
aparen como sumas y productos, se llama grado de la ecuación a la
mayor potencia de una derivada de y que aparece.
En este capítulo estamos interesados exclusivamente en ecuacio-
nes diferenciales resueltas (que no significa que ya estén resueltas,
sino que tienen una estructura concreta):
y′ = f ( x, y).
DEFINICIÓN 25. Un problema de condición inicial es una ecuación di-
ferencial junto con una condición inicial de la forma y( x0 ) = y0 , don-
de x0 , y0 ∈ R.
DEFINICIÓN 26. La solución general de una ecuación diferencial E es
una familia de funciones f ( x, c) donde c es una (o varias) constantes
tal que:
100 6. ECUACIONES DIFERENCIALES

• Cualquier solución de E tiene la forma f ( x, c) para algún c.

• Cualquier expresión f ( x, c) es una solución de E,
salvo para quizás un número finito de valores de c.
Lo primero que hay que tener en cuenta es que, si ya el problema
de calcular una primitiva de una función es complejo, integrar una
ecuación diferencial es, por lo general, imposible. Es decir, calcular la
función que cumple una determinada ecuación diferencial y su con-
dición inicial (o calcular la solución general) es un problema que no se
desea resolver en general. Lo que se desea es (sobre todo en ingeniería),
calcular una solución aproximada y conocer una buena cota del error
cometido al utilizar dicha solución en lugar de la “real”.

3. Discretización
En todo este capítulo, se supondrá que se tiene una función f ( x, y)
de dos variables, definida en una región x ∈ [ x0 , xn ], y ∈ [ a, b], que
cumple la siguiente condición:
DEFINICIÓN 27. Se dice que f ( x, y) definida en un conjunto X ∈ R2
cumple la condición de Lipschitz si existe una constante K > 0 tal que
| f ( x1 ) − f ( x2 )| ≤ K | x1 − x2 |
para cualesquiera x1 , x2 ∈ X, donde | | denota el módulo de un vec-
tor.
La condición de Lipschitz es una forma de continuidad un poco
fuerte (es decir, es más fácil ser continua que Lipschitz). Lo impor-
tante de esta condición es que asegura la existencia y unicidad de so-
lución de las ecuaciones diferenciales “normales”. Sea X un conjunto
de la forma [ x0 , xn ] × [ a, b] (una banda vertical) y f ( x, y) : X → R
una función de Lipschitz en X:
TEOREMA 9 (de Cauchy-Kovalevsky). Si f cumple las condiciones
anteriores, cualquier ecuación diferencial y′ = f ( x, y) con una condición
inicial y( x0 ) = y0 con y0 ∈ ( a, b) tiene una única solución y = y( x )
definida en [ x0 , x0 + t] para cierto t ∈ R mayor que 0.
No es difícil que una función cumpla la condición de Lipschitz:
basta con que sea, por ejemplo, diferenciable con continuidad en to-
dos los puntos y tenga derivada acotada. De hecho, todos los poli-
nomios y todas las funciones derivables que se han utilizado en este
curso son de Lipschitz en conjuntos de la forma [ x0 , xn ] × [ a, b], sal-
vo que tengan una discontinuidad o algún √ punto en el que no sean
derivables. Por ejemplo, la función f ( x ) = x no es de Lipschitz si el
 3. DISCRETIZACIÓN 101

intervalo [ x0 , xn ] contiene al 0 (porque el crecimiento de f cerca de 0

es “vertical”).
EJEMPLO 13 (Bourbaki “Functions
√ of a Real Variable”, Ch. 4, §1).
′
La ecuación diferencial y = 2 |y| con la condición inicial y(0) = 0
tiene infinitas soluciones. Por ejemplo, cualquier función definida así:
(1) y(t) = 0 para cualquier intervalo (−b, a),
(2) y(t) = −(t + β)2 para t ≤ −b,
(3) y(t) = (t − a)2 para t ≥ a
es una solución
√ de dicha ecuación diferencial. Esto se debe a que la
función |y| no es de Lipschitz cerca de y = 0. (Compruébense ambas
afirmaciones: que√dichas funciones resuelven la ecuación diferencial
y que la función |y| no es de Lipschitz).
Es decir, cualquier problema de condición inicial “normal” tiene
una única solución. Lo difícil es calcularla con exactitud.
¿Y cómo puede calcularse de manera aproximada?

3.1. La derivada como flecha. Es costumbre pensar en la deriva-

da de una función como la pendiente de la gráfica de dicha función en
el punto correspondiente. Es más útil pensar que es la coordenada Y
del vector velocidad de la gráfica.
Cuando se representa gráficamente una función, lo que se debe
pensar es que se está trazando una curva con velocidad horizontal
constante (el eje OX tiene la misma escala en todas partes, “se va de
izquierda a derecha con velocidad uniforme”). De esta manera, una
gráfica de una función f ( x ) es una curva ( x, f ( x )). El vector tangente
a esta curva es (derivando respecto del parámetro, que es x) el vec-
tor (1, f ′ ( x )): la derivada de f indica la dirección vertical del vector
tangente a la gráfica si la velocidad horizontal es constante.
De esta manera, se puede entender una ecuación diferencial de
la forma y′ = f ( x, y) como una ecuación que dice “se busca la curva
( x, y( x )) tal que el vector velocidad en cada punto es (1, f ( x, y)).” Así
que puede dibujarse el conjunto de vectores “velocidad” en el plano
( x, y), que están dados por (1, f ( x, y)) (la función f es conocida) y
se trata de dibujar una curva que tenga como vector velocidad ese en
cada punto y (la condición inicial) que pase por un punto concreto
(y( x0 ) = y0 ). Como se ve en la Figura 1, si se visualiza así, es sencillo
hacerse una idea de cómo será la curva que se busca.
Si se tuvieran las flechas (los vectores (1, f ( x, y))) dibujados en
un plano, trazar una curva tangente a ellos no es necesariamente muy
difícil. De hecho, si lo que se quiere es una aproximación, la idea más
102 6. ECUACIONES DIFERENCIALES

intuitiva es “sustituir la curva por segmentos que van en la dirección

de las flechas”. Si el segmento se toma con longitud muy pequeña,
uno piensa que se aproximará a la solución.

FIGURA 1. Flechas que representan vectores de una

ecuación diferencial del tipo y′ = f ( x, y). Obsérvese
que todos los vectores tienen misma magnitud en hori-
zontal.

Esta es precisamente la idea de Euler.

Ante un plano “lleno de flechas” que indican los vectores veloci-
dad en cada punto, la idea más simple para trazar una curva que
• Pase por un punto especificado (condición inicial y( x0 ) =
y0 )
• Tenga el vector velocidad que indica f ( x, y) en cada punto
La idea más sencilla es discretizar las coordenadas x. Partiendo de
x0 , se supone que el eje OX está cuantizado en intervalos de anchura
h —de momento una constante fija— “lo suficientemente pequeña”.
Con esta hipótesis, en lugar de trazar una curva suave, se aproxima la
solución dando saltos de anchura h en la variable x.
Como se exige que la solución pase por un punto concreto y( x0 ) =
y0 , no hay más que:
• Trazar el punto ( x0 , y0 ) (la condición inicial).
• Computar f ( x0 , y0 ), que es el valor vertical del vector veloci-
dad de la curva en cuestión en el punto inicial.
• Puesto que las coordenadas x están cuantizadas, el siguiente
punto de la curva tendrá coordenada x igual a x0 + h.
 4. [Contenido no esencial] ERRORES: RUNCAMIENTO Y REDONDEO 103

• Puesto que la velocidad en ( x0 , y0 ) es (1, f ( x0 , y0 )), la apro-

ximación más simple a cuánto se desplaza la curva en un inter-
valo de anchura h en la x es (h, h f ( x0 , y0 )) (tomar la x como el
tiempo y desplazarse justamente ese tiempo en línea recta).
Llámense x1 = x0 + h e y1 = y0 + h f ( x0 , y0 ).
• Trácese el segmento que une ( x0 , y0 ) con ( x1 , y1 ): este es el
primer “tramo aproximado” de la curva buscada.
• En este momento se tiene la misma situación que al principio
pero con ( x1 , y1 ) en lugar de ( x0 , y0 ). Repítase el proceso con
( x1 , y1 ).
Y así hasta que se desee. Lo que se acaba de describir se conoce co-
mo el algoritmo de Euler para integrar numéricamente una ecuación
diferencial.

4. [Contenido no esencial] Errores: runcamiento y redondeo

Obviamente, la solución de una ecuación diferencial y′ = f ( x, y)
no será siempre (ni casi nunca) una curva formada por segmentos
rectilíneos: al aproximar la solución y( x ) por los segmentos dados
por el método de Euler, se está cometiendo un error que puede ana-
lizarse con la fórmula de Taylor. Si se supone que f ( x, y) es suficien-
temente diferenciable, entonces y( x ) también lo será y se tendrá que

h2 ′′
y( x0 + h) = y( x0 ) + hy′ ( x0 ) + y (θ )
2
para cierto θ ∈ [ x0 , x0 + h]. Como lo que se está haciendo es olvidarse
del último sumando, el error que se produce es exactamente ese suman-
do: un término del tipo O(h2 ), de segundo orden al hacer la primera
iteración. Este error, que aparece al aproximar una función por su de-
sarrollo limitado de Taylor de un grado finito, se denomina error de
truncamiento (pues se está truncando la expresión exacta de la fun-
ción solución).
Por lo general, se supone siempre que el intervalo en que se inte-
gra la ecuación diferencial es del orden de magnitud de h−1 . De he-
cho, si el intervalo de cálculo es [ a, b], se fija un número de intervalos n
y se toma x0 = a, xn = b y h = (b − a)/n. Por tanto, realizar el cálculo
desde x0 = a hasta xn = b requiere n = 1/h iteraciones, de manera
que el error de truncamiento global es del orden de h−1O(h2 ) ≃ O(h):
tiene magnitud similar a la anchura intervalo. Aunque, como se ve-
rá, esto hay que entenderlo bien, pues por ejemplo (b − a)10 h es del
orden de h, pero si b ≫ a, entonces es mucho mayor que h. En estos
104 6. ECUACIONES DIFERENCIALES

problemas en que se acumulan errores, es importante tener en cuenta

que la acumulación puede hacer aparecer constantes muy grandes.
Pero el truncamiento no es la única parte de la que surgen erro-
res. Como ya se ha insistido, calcular en coma flotante (que es lo
que se hace en este curso y en todos los ordenadores, salvo raras ex-
cepciones) lleva implícito el redondeo de todos los datos a una pre-
cisión específica. En concreto, si se utiliza coma flotante de doble
precisión como especifica el estándar IEEE-754, se considera que el
número más pequeño significativo (el épsilon de dicho formato) es
ϵ = 2−52 ≃ 2.2204 × 10−16 . Es decir, lo más que se puede suponer es
que el error de redondeo es menor que ϵ en cada operación. Si se supo-
ne que “operación” significa “iteración del algoritmo”, se tiene que
cada paso de xi a xi+1 genera un error de magnitud ϵ. Como hay que
realizar n = 1/h pasos desde x0 hasta xn , el error de redondeo global es
del orden de ϵ/h. Por tanto, como ϵ es un dato fijo (no se puede mejo-
rar en las circunstancias de trabajo), cuanto más pequeño sea el intervalo
entre xi y xi+1 , mayor será el error de redondeo.
Así pues, la suma del error de truncamiento y el de redondeo es
de la forma (tomando infinitésimos equivalentes)
ϵ
E(h) ≃ + h,
h
que crece cuando h se hace grande y cuando h se hace pequeño.√El
mínimo de E(h) se alcanza (esta es una fácil cuenta) cuando √ h ≃ ϵ:
no tiene sentido utilizar intervalos de anchura menor que ϵ porque
el error de redondeo se magnificará. Por tanto: tomar intervalos lo más
pequeños posibles puede conducir a errores de gran magnitud.
En el caso de doble precisión, la anchura más pequeña sensata
es del orden de 10−8 : anchuras menores son inútiles, además de la
sobrecarga de operaciones que exigen al ordenador. Si en algún mo-
mento parece necesario elegir una anchura menor, lo más probable
es que lo que se requiera es elegir un algoritmo mejor.
Esta consideración sobre el error de redondeo y de truncamiento
es independiente del método: cualquier algoritmo discreto incurrirá
en un error de truncamiento inherente al algoritmo y en un error de
redondeo. Lo que se ha de hacer es conocer cómo se comporta cada
algoritmo en ambos casos.

4.1. Anchura variable. Hasta ahora se ha supuesto que el inter-

valo [ a, b] (o, lo que es lo mismo, [ x0 , xn ]) se divide en subintervalos
de la misma anchura h. Esto no es necesario en general. Pero simpli-
fica la exposición lo suficiente como para que siempre supongamos
 6. EL MÉTODO DE EULER 105

que los intervalos son de anchura constante, salvo que se indique ex-
plícitamente lo contrario (y esto es posible que no ocurra).

5. [Contenido no esencial] Solución de EDOs e integral

Supóngase que la ecuación diferencial y′ = f ( x, y) es tal que la
función f depende solo de la variable x. En este caso, la ecuación se
escribirá
y′ = f ( x )
y significa y es la función cuya derivada es f ( x ), es decir, es el problema
de calcular una primitiva. Este problema ya se ha estudiado en su
capítulo, pero está intrínsecamente ligado a casi todos los algoritmos
de integración numérica y por ello es relevante.
Cuando la función f ( x, y) depende también de y, la relación es
más difícil, pero también se sabe (por el Teorema Fundamental del
Cálculo) que, si y( x ) es la solución de la ecuación y′ = f ( x, y) con
condición inicial y( x0 ) = y0 , entonces, integrando ambos miembros,
queda
∫ x
y( x ) = f (t, y(t)) dt + y0 ,
x0
que no es el cálculo de una primitiva pero se le parece bastante. Lo que
ocurre en este caso es que los puntos en los que hay que evaluar
f ( x, y) no se conocen, pues justamente dependen de la gráfica de la
soclución (mientras que al calcular una primitiva, el valor de f ( x, y)
no depende de y). Lo que se hace es intentar acertar lo más posible con la
información local que se tiene.

6. El método de Euler: integral usando el extremo izquierdo

Si se explicita el algoritmo para integrar una ecuación diferencial
utilizando el método de Euler como en el Algoritmo 10, se observa
que se aproxima cada valor siguiente yi+1 como el valor anterior al
que se le suma el valor de f en el punto anterior multiplicado por la
anchura del intervalo. Es decir, se está aproximando
∫ x i +1
f (t, y(t)) dt ≃ ( xi+1 − xi ) f ( xi , yi ) = h f ( xi , yi ).
xi

Si la función f ( x, y) no dependiera de y, se estaría haciendo la aproxi-

mación
∫ b
f (t) dt = (b − a) f ( a),
a
106 6. ECUACIONES DIFERENCIALES

que se puede denominar (por falta de otro nombre mejor) la regla del
extremo izquierdo: se toma como área bajo f ( x ) la anchura del intervalo
por el valor en el extremo izquierdo de este.

Algoritmo 10 (Algortimo de Euler con aritmética exacta)

Input: Una función f ( x, y), una condición inicial ( x0 , y0 ), un inter-
valo [ x0 , xn ] y un paso h = ( xn − x0 )/n
Output: una colección de puntos y0 , y1 , . . . , yn (que aproximan los
valores de la solución de y′ = f ( x, y) en la malla x0 , . . . , xn )
⋆INICIO
i ← 0, n ← ( xn − x0 )/h
while i ≤ n do
y i +1 ← y i + h f ( x i , y i )
i ← i + 1, xi ← xi−1 + h
end while
return (y0 , . . . , yn )

Podría intentarse (se deja de momento como ejercicio) plantearse

el problema “y si se quisiera utilizar el extremo derecho del intervalo,
¿qué habría que hacer?”: esta pregunta da lugar al método implícito
más simple, que todavía no estamos en condiciones de explicar.

7. [Contenido no esencial] Euler modificado: “punto medio”

En lugar de utilizar el extremo izquierdo del intervalo [ xi , xi+1 ]
para “integrar” y calcular el siguiente yi+1 , sería preferible intentar
utilizar la regla del punto medio como primera idea. El problema es
que con los datos que se tienen, no se conoce cuál sea dicho punto
medio (por la dependencia de y que tiene f ). Aun así, se puede tratar
de aproximar de la siguiente manera:
• Se utilizará un punto cercano a Pi = ( xi , yi ) y cuya coorde-
nada x esté mitad de la distancia entre xi y xi+1 .
• Como no se conoce nada mejor, se hace la primera aproxi-
mación como en el algoritmo de Euler y se toma como ex-
tremo derecho Qi = ( xi+1 , yi + h f ( xi , yi )).
• Se calcula el punto medio del segmento Pi Qi , que es ( xi +
h/2, yi + h/2 f ( xi , yi )). Llámese k a su coordenada y.
• Se utiliza el valor de f en este punto para calcular el siguiente
yi+1 : yi+1 = yi + h f ( xi + h/2, k ).
 7. [Contenido no esencial] EULER MODIFICADO: “PUNTO MEDIO” 107

Si f ( x, y) no depende de y, es sencillo comprobar que lo que se

está haciendo con los pasos descritos es la aproximación
∫ b
a+b
f ( x ) dx ≃ (b − a) f ( ),
a 2
que es la regla del punto medio de integración numérica.
Este método se denomina método de Euler modificado y tiene un
orden de precisión superior al de Euler; es decir, es de orden 2, lo que
significa que el error acumulado en xn es O(h2 ). Se describe en detalle
en el Algoritmo 11.

yi + k1
z2 + k 2

z2
y i +1

yi

xi xi + h x i +1
2

FIGURA 2. Ilustración de Algoritmo de Euler modifica-

do. En lugar de sumar el vector (h, h f ( xi , yi )) (azul),
se suma en el punto ( xi , yi ) el vector verde, que corres-
ponde al vector tangente dado en el punto medio (en
rojo).

Como muestra la Figura 2, el algoritmo de Euler modificado con-

siste en utilizar, en lugar del vector tangente a la curva descrito por
(1, f ( xi , yi )), “llevar a ( xi , yi )” el vector tangente en el punto medio
del segmento que une ( xi , yi ) con el punto que correspondería a utili-
zar el algoritmo de Euler. Desde luego, se sigue cometiendo un error,
pero al tomar la información un poco más lejos, se aproxima algo mejor
el comportamiento de la solución real.
La siguiente idea es utilizar una “media” de vectores. En concreto,
la media entre los dos vectores, el del “inicio del paso” y el del “final
del paso de Euler.”
108 6. ECUACIONES DIFERENCIALES

Algoritmo 11 (Algortimo de Euler Modificado con aritmética exacta)

Input: Una función f ( x, y), una condición inicial ( x0 , y0 ), un inter-
valo [ x0 , xn ] y un paso h = ( xn − x0 )/n
Output: una colección de puntos y0 , y1 , . . . , yn (que aproximan los
valores de la solución de y′ = f ( x, y) en la malla x0 , . . . , xn )
⋆INICIO
i←0
while i ≤ n do
k1 ← f ( xi , yi )
z2 ← yi + 2h k1
k2 ← f ( xi + 2h , z2 )
yi+1 ← yi + hk2
i ← i+1
end while
return (y0 , . . . , yn )

8. Método de Heun (regla del trapecio)

En lugar de utilizar el punto medio del segmento que une el inicio
y el final del método de Euler puede hacerse una operación parecida
con los vectores: utilizar la media entre el vector al inicio del paso y el
vector al final del paso de Euler. Este es el algoritmo de Euler mejorado
ó algoritmo de Heun. Se describe en detalle en el Algoritmo 12.
Es decir, cada paso consiste en las siguientes operaciones:
• Se calcula el valor k1 = f ( xi , yi ).
• Se calcular el valor z2 = y j + hk1 . Esta es la coordenada yi+1
en el algoritmo de Euler.
• Se calcula el valor k2 = f ( xi+1 , z2 ). Esta sería la pendiente
en el punto ( xi+1 , yi+1 ) en el algoritmo de Euler.
1 2
• Se calcula la media de k1 y k2 : k +2 k y se usa este valor como
“paso”, es decir: yi+1 = yi + 2h (k1 + k2 ).
El método de Heun se describe gráficamente en la Figura 3.

9. El modelo Predador-Presa (o de Lotka-Volterra)

El primer sistema de ecuaciones diferenciales que modela un sis-
tema biológico complejo es el de “Lotka-Volterra”: se utiliza para des-
cribir la evolución de un entorno en el que dos especies conviven; una
de ellas actúa como depredador y otra como presa (el ejemplo más
habitual que se da es el de zorros y conejos).
 9. EL MODELO PREDADOR-PRESA (O DE LOTKA-VOLTERRA) 109

yi + hk1

y i +1

yi + hk2

yi
xi x i +1 x i +1 + h

FIGURA 3. Ilustración de Algoritmo de Heun. En el pun-

to ( xi , yi ) se añade la media (en verde) de los vectores
correspondientes al punto ( xi , yi ) (azul) y al siguiente
de Euler (en rojo).

Algoritmo 12 (Algortimo de Heun con aritmética exacta)

Input: Una función f ( x, y), una condición inicial ( x0 , y0 ), un inter-
valo [ x0 , xn ] y un paso h = ( xn − x0 )/n
Output: una colección de puntos y0 , y1 , . . . , yn (que aproximan los
valores de la solución de y′ = f ( x, y) en la malla x0 , . . . , xn )
⋆INICIO
i ← 0, n ← ( xn − x0 )/h
while i ≤ n do
m e ← f ( xi , yi )
ye ← yi + hk1
m f ← f ( xi + h, z2 )
yi+1 ← yi + 2h (me + m f )
i ← i + 1, xi ← xi−1 + h
end while
return (y0 , . . . , yn )

Supongamos que x (t) denota la población de la especie “presa”

en el momento t y que y(t) denota la población de la especie “preda-
dor”. La ecuación diferencial de Lotka-Volterra que describe la evo-
lución conjunta de las dos especies se basa en las siguientes suposi-
ciones (simplistas, claro está):
110 6. ECUACIONES DIFERENCIALES

(1) La población de presas crece, sin mirar la interacción con los

predadores, en proporción a su tamaño, pues tiene alimento
suficiente (pongamos, yerba). En realidad, las muertes están
incluidas en este apartado, pues basta reducir el incremento
para incluirlas.
(2) Las presas, además, mueren como consecuencia de ser víc-
tima de un depredador. Esto ocurre con una probabilidad
constante.
(3) La población de depredadores muere, en ausencia de inter-
acción con las presas, en proporción a su tamaño.
(4) Los predadores solo se multiplican en proporción a las pre-
sas que comen.
Estas cuatro reglas elementales (insistimos, muy simples), se tras-
ladan a ecuaciones como se explica a continuación. Pese a lo elemen-
tal de la descripción, el modelo, como se verá presenta propiedades
interesantes.
Del punto 1 se deduce que existe una constante α > 0 tal que
ẋ (t) = αx (t) + . . .
donde los puntos denotan otra expresión que habrá que descubrir. Si
solo se tuviera esta ecuación, habría un crecimiento de presas expo-
nencial (la ecuación ẋ (t) = αx (t) tiene como solución x (t) = Keαt ,
donde K es el valor inicial).
Del punto 3 se deduce que existe un número γ > 0 tal que
ẏ(t) = −γy(t) + . . .
igual que antes para las presas. Esto hace que, de por sí, los depreda-
dores decrezcan según una ley exponencial negativa.
Es fácil convencerse de que la probabilidad de que un depredador
se encuentre con una presa en algún lugar del terreno es proporcio-
nal al producto del número de predadores y de presas. Pero como
no todo encuentro mutuo termina en una caza, se modela la proba-
bilidad de que un predador coma a una presa como una constante β
por dicho producto: βx (t)y(t) (con β > 0). Por otro lado, el hecho de
que un depredador se alimente no garantiza que se reproduzca; hare-
mos que esto ocurra con un factor δ. Así, los puntos suspensivos que
hemos dejado arriba han de sustituirse por − βx (t)y(t) (en la parte
de las presas) y por δx (t)y(t) (en la parte de los depredadores). Se
obtiene el sistema de ecuaciones diferenciales de Lotka-Volterra:
ẋ (t) = αx (t) − βx (t)y(t)
(18)
ẏ(t) = −γy(t) + δx (t)y(t)
 9. EL MODELO PREDADOR-PRESA (O DE LOTKA-VOLTERRA) 111

9.1. Cómo se aplica Heun al sistema Predador-Presa. En lugar

de tratar de implementar el algoritmo de Euler y luego el de Heun
para un sistema de ecuaciones, describimos directamente el segun-
do (así que el primero queda incluido en él). Como ya se explicó, el
algoritmo de Euler es casi siempre inadecuado.
En el momento en que aparece más de una variable, los cálculos
se multiplican pero la filosofía es siempre la misma: discretizar la va-
riable independiente e ir paso a paso sin precipitarse.
El Problema de Condición Inicial es:
{
ẋ = 0.8x − 0.4xy x (0) = 18
ẏ = −2y + 0.2xy y(0) = 3
(los valores iniciales en un sistema Predador-Presa siempre se conci-
ben como densidades de población, no como números absolutos, pues si
fuera esto llegaríamos a poblaciones de “medio zorro” o “doscientos
conejos y cuarto”).
La variable independiente es el tiempo, que denotamos con una
t. Fijemos una discretización, por ejemplo h = 0.25 y calculemos dos
pasos (más sería demasiado pesado).
La única diferencia con el caso de una variable es que, en cada pa-
so, hay que hacer tantos cálculos como variables dependientes. Dos
para este ejemplo.
La función que da la derivada de x con respecto a t es f (t, x, y) =
0.8x − 0.4xy. Para la variable y, la función correspondiente es g(t, x, y) =
−2y + 0.2xy. Obsérvese —esto es importante— que tanto f como g po-
drían depender de t, aunque en este caso no ocurra así.
(1) Primer paso, t0 = 0. Tenemos que x̃0 = x (0) = 18 y que
ỹ0 = y(0) = 3.
(a) La pendiente de Euler de la x en el punto inicial es:
me,x = f (t0 , x̃0 , ỹ0 ) = −7.2.
(b) La pendiente de Euler de la y en el punto inicial es:
me,y = g(t0 , x̃0 , ỹ0 ) = 4.8.
(c) Calculamos la predicción que daría el algoritmo de Eu-
ler para la x:
xe = x̃0 + h · me,x = 18 + 0.25 · (−7.2) = 16.2.
(d) Calculamos la predicción que daría el algoritmo de Eu-
ler para la y:
ye = ỹ0 + h · me,y = 3 + 0.25 · 4.8 = 4.2.
112 6. ECUACIONES DIFERENCIALES

(e) Calculamos la pendiente de la x en el punto final:

m f ,x = f (t0 , x̃e , ỹe ) = −14.256.
(f) Calculamos la pendiente de la y en el punto final:
m f ,y = g(t0 , x̃e , ỹe ) = 5.208.
(g) Las pendientes medias son, por tanto:
m x = (−7.2 − 14.256)/2 = −10.728,
my = (4.8 + 5.208)/2 = 5.004.
(h) Así que, por fin:
x1 = x0 + h · m x = 15.318
y1 = y0 + h · my = 4.251.
(2) Segundo paso, t1 = 0.25. Tenemos que x̃1 = 15.318 y que
ỹ1 = 4.251.
(a) Pendientes de Euler (las dos):
me,x = f (t1 , x1 , y1 ) = −13.79233
me,y = g(t1 , x1 , y1 ) = 4.5214.
(b) Predicciones de Euler:
xe = x1 + h · me,x = 11.87
ye = y1 + h · me,y = 5.381.
(c) Pendientes en el punto final:
m f ,x = f (t1 , xe , ye ) = −16.055
m f ,y = g(t1 , xe , ye ) = 2.013.
(d) Pendientes medias:
m x = (me,x + m f ,x )/2 = −14.923
my = (me,y + m f ,y )/2 = 3.267.
(e) Valores de x̃2 , ỹ2 :
x̃2 = x̃1 + h · m x = 11.587
ỹ2 = ỹ1 + hṁy = 5.068.
Los valores exactos son x (0.5) = 11.7 e y(0.5) = 5.131. Se aprecia una
muy buena aproximación para lo grande que es el paso (0.25).
Pero lo más claro es lo laborioso del cálculo. Hoy día estas cuentas
se realizan por ordenador. Durante el programa Apollo, la mayoría
de estas operaciones se llevaban a cabo a mano. Quienes lo hacían
 9. EL MODELO PREDADOR-PRESA (O DE LOTKA-VOLTERRA) 113

son los héroes no reconocidos de la Carrera Espacial (aunque ya se

ha hecho alguna película al respecto).
9.2. Representación de una solución. En la Figura 4, se ha repre-
sentado un ejemplo de solución (usando el algoritmo de Heun) del
sistema de Lotka-Volterra de la sección anterior. Obsérvese (esto es
quizás lo más importante de este sistema) que ambas poblaciones tie-
nen un comportamiento creciente y decreciente con un desfase entre
ellas.
Con algo de esfuerzo se puede verificar que los puntos críticos de
la función x (t) se alcanzan cuando la función y(t) vale γ/δ y que los
puntos críticos de y(t) se alcanzan cuando la función x (t) vale α/β
(¿cómo se puede comprobar esto? es fácil pero requiere una explica-
ción).
Sin embargo, se sabe (desde el punto de vista teórico) que el siste-
ma de Lotka-Volterra es periódico. Las soluciones que se han dibujado
no lo son (si uno calculara las soluciones para tiempos mucho más
largos, vería cómo las funciones se vuelven cada vez más exagerada-
mente “puntiagudas”).

conejos
zorros

0 5 10 15 20

FIGURA 4. Ejemplo de sistema predador-presa. En este

sistemas se supone siempre que los valores son “densi-
dades de población”.