Tema 4

Interpolación

4.1. Introducción
El problema matemático de la interpolación es el siguiente: Dada un conjunto de n pares de
valores (xk , yk ), encontrar una función f (x) que cumpla f (xk ) = yk , k = 1, n. Existen diversos
métodos para encontrar dicha función. Los más conocidos son los métodos que interpolan f (x)
mediante un polinomio o una funcion racional (cociente de dos polinomios). En los últimos años.
los métodos que utilizan funciones definidas a tramos (splines) han cobrado gran popularidad,
justificada por su capacidad de reproducir formas complicadas y altamente variables. De hecho
existe una activa investigación es splines.
Vamos a estudiar en primer lugar la interpolación polinómica. En este caso, dados n puntos,
el grado más bajo del polinomio que pasa por los n puntos es n − 1 (salvo que dos o más puntos
pertenezcan a un polinomio de grado más bajo). El Teorema Fundamental del Algebra garantiza
que este polinomio de grado n − 1 es único. Existen diversos métodos de determinar el polinomio
interpolador, que dan por supuesto el mismo resultado. Unos métodos son más conveniente que
otros, dependiendo del número de puntos en los que se desee conocer el polinomio interpolador.

4.2. Método de interpolación de Lagrange

Este método consiste en construir el polinomio interpolador de grado n que pasa por n + 1
puntos (xi , yi ) de la forma
n
Pn (x) = ∑ Li (x)yi
i=0
donde las funciones Li (x) cumplen Li (xk ) = 0 si i "= k y Li (xi ) = 1. Esta propiedad garantiza
Pn (xk ) = yk . Las funciones Li (x) se construyen como:

(x − x1 )(x − x2 ) · · · (x − xi−1 )(x − xi+1 ) · · · (x − xn )

Li (x) =
(xi − x1 )(xi − x2 ) · · · (xi − xi−1 )(xi − xi+1 ) · · · (xi − xn )

58
CAPÍTULO 5. INTERPOLACIÓN 59

De su forma explícita, vemos que Li (xi ) = 1 y que Li (xk ) = 0 para i "= k. Podemos escribir el
polinomio interpolador de una función f (x) en los puntos x0 , x1 , . . . , xn de forma compacta como:
n n (x − x j )
Pn (x) = ∑ f (xi ) ∏ (xi − x j )
i=0 j=0
j "= i

El método de Lagrange permite construir fácilmente de forma explícita, el polinomio inter-

polador. Dos de las interpolaciones más utilizados son las interpolaciones lineal y parabólica.
Los correspondientes polinomios interpoladores son:
Lineal
x − x0 x − x1
P(x) = f1 + f0
x1 − x0 x0 − x1
Parabólica

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

P2 (x) = f0 + f1 + f2
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )
Sin embargo, cuando deseamos evaluar el polinomio interpolador en uno o varios puntos, el
método de Lagrange no es conveniente ya que no da los coeficiente de cada grado, lo cual es
necesario para aplicar el algoritmo de Horner. Sin embargo el método de Lagrange es fácil de
entender, lo que justifica su valor pedagógico, y es de utilidad a la hora de elaborar un cierto
número de demostraciones, como por ejemplo la fórmula del error cometido por el polinomio
interpolador, que pasamos a demostrar a continuación.
Un aspecto importante es la determinación del error del polinomio interpolador Pn (x) con
respecto a la función interpolada f (x), en el intervalo [x0 , xn ]. Sea xs un punto arbitrario de este
intervalo, en el que deseamos estimar el error. Construimos la función
n (x − x j )
g(x) = f (x) − Pn (x) − ∏ ( f (xs ) − Pn (xs ))
(xs − x j )
j=0

que, por construcción, se anula en x0 , x1 , . . . xn y xs . Como g(x) tiene n + 2 ceros y por el teorema
de Rolle la derivada g(1) (x) tiene n + 1 ceros, la derivada g(2) (x) tiene n ceros y así sucesivamente
hasta la derivada g(n+1) que tiene al menos un cero. Sea α este cero. Derivando n + 1 veces la
anterior igualdad y sustituyendo x = α obtenemos (la derivada n + 1 de Pn (x) es nula y la del
polinomio de orden n + 1 del numerador del último término es (n + 1)!)
(n + 1)!
g(n+1) (α) = 0 = f (n+1) (α) − ( f (xs ) − Pn (xs )
∏ n
(xs − x j )
j=0
con lo cual, concluimos que:
n
f (n+1) (α)
f (xs ) = Pn (xs ) + ∏ (xs − x j )
(n + 1)!
j=0
CAPÍTULO 5. INTERPOLACIÓN 60

El punto α es desconocido, pero si podemos acotar la función f (n+1) (x) en el intervalo [x0 , xn ] por
el valor f (n+1) (xm ), donde xm es el punto donde f (n+1) (x) alcanza el valor máximo, encontramos
la cota de error en el intervalo [x0 , xn ]:
! !
! n (n+1)
!
!∏ (x − x ) f (x )!
! j=0 j m !
| f (x) − Pn (x)| ≤
(n + 1)!
que es válida, independientemente del método empleado para obtener el polinomio interpolador.

4.3. Algoritmo de Neville

A veces no es necesario obtener la forma exolícita del polinomio interpolador y basta con
obtener su valor numérico en un punto dado. Además en este caso nos gustaría el poder aumen-
tar el orden del polinomio interpolador a voluntad y parar cuando el error sea suficientemente
pequeño. Para estos propósitos el algoritmo de Neville está especialmente indicado. El algoritmo
de Neville está inspirado en la fórmula de interpolación lineal. La expresión
x − x0 x − xk
P(x) = fk − f0
xk − x0 xk − x0
pasa por los puntos (x0 , f0 ) y (xk , fk ). Si, en lugar de fk , ponemos en la anterior expresión un
(k−1)
polinomio interpolador P12···k (x) que pasa por una serie de puntos xr , r = 1, 2, . . . , k, y en lugar
de f0 , un polinomio interpolador P01···k−1 (x) que pasa por (xr , fr ), r = 0, 1, . . . , k − 1, obtenemos
x − x0 (k−1) x − xk (k−1)
P(x) = P12···k (x) − P (x)
xk − x0 xk − x0 01···k−1
Es obvio de la expresión anterior que, para los puntos comunes de interpolación xr , r = 1, 2, . . . , k-
(k−1) (k−1)
1 de P0 (x) y Pk (x), se cumple
xr − x0 xr − x1
P(xr ) = fr − fr = fr
x1 − x0 x1 − x0
(k−1) (k−1)
y que además, P(x0 ) = P01···k−1 (x0 ) = f0 y P(xk ) = P12···k (xk ) = f1 . Por lo tanto podemos cons-
truir un polinomio de orden k a partir de dos polinomios de orden k − 1, mediante la fórmula de
interpolación lineal. El algoritmo de Neville utiliza este hecho, conjuntamente con una adecuada
redisposición de los puntos de interpolación en una tabla, para calcular el valor del polinomio
interpolador en un punto x. La tabla que se construye, en el caso de cuatro puntos, es
x − x0 f0
(1)
P01
(2)
x − x1 f1 P012
(1) (3)
P12 P0123
(2)
x − x2 f2 P123
(1)
P23
x − x3 f3
CAPÍTULO 5. INTERPOLACIÓN 61

donde cada columna se obtiene de la anterior mediante la fórmula de interpolación lineal, que
equivale a un producto cruzado entre la primera columna y la columna previa a la que calculamos.
Concretamente:
(1) 1
P01 (x) = [(x − x0 ) f1 − (x − x1 ) f0 ]
x1 − x0
(2) 1 (1) (1)
P012 (x) = [(x − x0 )P12 (x) − (x − x1 )P01 (x)]
x2 − x0
(3) 1 (2) (2)
P0123 (x) = [(x − x0 )P123 (x) − (x − x3 )P012 (x)]
x3 − x0
Vemos que cada término de la tabla se obtiene como la diferencia de los productos cruzados
de vecinos inferior y superior en la columna inmediatamente anterior con la diferencia entre la
abcisa de evaluación y la abcisa de interpolación en la primera columna correspondiente con el
otro vecino, siguiendo la diagonal, dividida (la diferencia de productos cruzados) por la diferen-
cia entre las diferencias de abcisas correspondientes superior e inferior de la primera columna.
Podemos escribir el término general.como sigue. Si denominamos y j = x − x j , tenemos:
1 " #
(i+1) (i) (i)
Ps···(s+i+1) = ys Ps+1···s+i+1 − ys+i+1 Ps···s+i
ys − ys+i+1
El algoritmo de Nevillel es particularmente indicado en el caso de que se desee evaluar el poli-
nomio interpolador en un único punto, o en un número muy reducido de puntos.

4.4. Método de Newton o de diferencias divididas

El método de Newton de diferencias divididas es otra forma de obtener el polinomio interpo-
lador. En este método el polinomio interpolador se escribe de la forma
Pn (x) = a0 + (x − x0 )a1 + (x − x0 )(x − x1 )a2 + · · · + (x − x0 )(x − x1 ) · · · (x − xn−1 )an
y el algoritmo proporciona una regla para obtener los coeficientes a0 , a1 , . . . , an . Imponiendo que
el polinomio interpolador pase por los puntos de interpolación obtenemos
Pn (x0 ) = a0 = f (x0 )
Pn (x1 ) = a0 + (x1 − x0 )a1 = f (x1 )
Pn (x2 ) = a0 + (x2 − x0 )a1 + (x2 − x0 )(x2 − x1 )a2 = f (x2 )
.. .. ..
. . .
Pn (xn ) = a0 + (xn − x0 )a1 + · · · + (xn − x0 ) · · · (xn − xn−1 )an = f (xn )
De estas ecuaciones, es obvio que a0 depende sólo de x0 , a1 de x0 y x1 y así sucesivamente.
Introducimos la nueva notación a0 ≡ f [x0 ], a1 ≡ f [x0 , x1 ], y así sucesivamentre, con f [x0 ] =
f (x0 ), como se ve de la primera ecuación. Restando las dos primeras ecuaciones obtenemos
f (x1 ) − f (x0 )
a1 ≡ f [x0 , x1 ] =
x1 − x0
CAPÍTULO 5. INTERPOLACIÓN 62

Restando la segunda y la tercera ecuación obtenemos:

f (x2 ) − f (x1 ) − f [x0 , x1 ](x2 − x1 ) f [x1 , x2 ] − f [x0 , x1 ]
a2 ≡ f [x0 , x1, x2 ] = =
(x2 − x0 )(x2 − x1 ) x2 − x0
Podemos proceder de igual modo para demostrar que
f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ]
an ≡ f [x0 , x1 , . . . , xn ] =
xn − x0
aunque la forma más cómoda es por inducción. Suponemos que la expresión vale para an−1 y
construimos el polinomio de grado n, Qn (x), definido por
% x − xn %
Qn (x) = Pn−1 (x) + [P (x) − Pn−1 (x)]
xn − x0 n−1
donde Pn−1% (x) es el polinomio interpolador en x , . . . , x . Está claro de la definición de Q (x)1
1 n n
que para x1 , x2 , . . . , xn−1 , Qn (xi ) = Pn−1
% (x ) = f (x ), ya que P% (x) − P
i i n−1 n−1 (x) se anula en estos
puntos. También en xn se cumple que Q(xn ) = Pn−1 (xn ) = f (xn ). En x0 se cumple Qn (x0 ) =
%
% (x ) − [P% (x ) − P
Pn−1 0 n−1 0 n−1 (x0 ) = Pn−1 (x0 ) = f (x0 ). Luego Qn (x) coincide con Pn (x), ya que el
polinomio interpolador es único, y su coeficiente en xn es por lo tanto an . Como hemos supuesto
que la fórmula de diferencias divididas es válida para Pn−1 (x) y Pn−1 % (x), entonces tenemos,

identificando potencias en xn en ambos lados, que que el coeficiente en xn de Qn (x) viene dado
por:
1 $ % % f [x1 , x2 , . . . , xn ] − f [x0 , x1 , . . . , xn−1 ]
an ≡ f [x0 , x1 , . . . , xn ] = an−1 − an−1 =
xn − x0 xn − x0
relación que es el origen del nombre de diferencias divididas para los coeficientes an .
Podemos por lo tanto escribir el polinomio interpolador como,
Pn (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + · · · f [x0 , x1 , . . . , xn ](x − x0 ) · · · (x − xn−1 )
y la fórmula de aproximación a f (x) con su término de error queda en la siguiente forma:
n
f (n+1) (xm )
f (x) = f [x0 ]+ f [x0 , x1 ](x−x0 )+· · ·+ f [x0 , x1 , . . . , xn ](x−x0 ) · · · (x−xn−1 )+ ∏ (x−x j )
(n + 1)!
j=0
El método de Newton permite obtener los coeficientes del polinomio interpolador fácilmente en
forma de tabla, que damos abajo para el caso de 4 puntos.
x0 f0
f [x0 , x1 ]
x1 f1 f [x0 , x1 , x2 ]
f [x1 , x2 ] f [x0 , x1 , x2 , x3 ]
x2 f2 f [x1 , x2 , x3 ]
f [x2 , x3 ]
x3 f3
1 Q (x) es el resultado de aplicar el algoritmo de Neville a Pn−1 (x) y Pn−1
% (x).
n
CAPÍTULO 5. INTERPOLACIÓN 63

El método de Newton es especialmente indicado en el caso de que deseemos realizar muchas

evaluaciones del polinomio interpolador, ya que da el polinomio preparado para ser evaluado por
el algoritmo de Horner. Otro aspecto particularmente conveniente es que, si deseamos aumentar
el orden del polinomio interpolador, los coeficientes ak ya calculados permanecen inalterados, es
decir, no destruímos el trabajo ya realizado cuando deseamos aumentar el orden del polinomio
interpolador. Se dice en este caso que los coeficientes ak tienen la propiedad de permanencia.
En el caso de puntos igualmente espaciados, el polinomio de Newton toma una forma espe-
cialmente conveniente. Supongamos que tenemos una red de puntos espaciados un paso h, de
forma que xn = x0 + nh. Tenemos que
f1 − f0
f [x0 , x1 ] =
h
f [x1 , x2 ] − f [x0 , x1 ] f2 − 2 f1 + f0
f [x0 , x1 , x2 ] = =
2h 2h2
Si introducimos la notación de diferencias finitas & f0 = f1 − f0 , &2 f0 = &(& f0 ) = f2 − f1 −
( f1 − f0 ) = f2 − 2 f1 + f0 llegamos fácilmente por inducción al resultado
&n f0
f [x0 , x1 , . . . , xn ] =
n!hn
ya que
&n−1 f1 &n−1 f0
f [x1 , . . . , xn ] − f [x0 , . . . , xn−1 ] (n−1)!hn−1
− (n−1)!hn−1 &(&n−1 f0 ) &n f0
f [x0 , x1 , . . . , xn ] = = = =
xn − x0 nh n!hn n!hn
El polinomio interpolador adquiere una forma particularmente simple en el caso de puntos
igualmente espaciados. Si denotamos un punto arbitrario x, comprendido entre x0 y xn , como
x0 + sh, tenemos que el término j del polinomio interpolador se puede expresar como:
& '
& j f0 s!h j & j f0 s
(x − x0 ) · · · (x − x j−1 ) f [x0 , · · · , x j ] = sh · · · (s − j + 1)h j
= j
= & j f0
j!h j!(s − j)! h j

Con esta extensión de los números combinatorios a números reales, podemos expresar el polino-
mio interpolador con su término de error en la siguiente forma extraordinariamente compacta:
n & ' & '
s s
Pn (x) = Pn (x0 + sh) = ∑ j
& f0 + hn+1 f (n+1) (ξ)
j=0
j n + 1

5.4.1. Relación de las derivadas con las diferencias

Las derivadas de f (x) y las diferencias finitas están íntimamente relacionadas. Para demos-
trarlo, vamos primero a expresar la función f (x) en término de diferencias divididas. La expre-
sión
f (x0 ) − f (x)
f [x, x0 ] =
x0 − x
CAPÍTULO 5. INTERPOLACIÓN 64

podemos reescribirla como

f (x) = f [x0 ] + (x − x0 ) f [x, x0 ]
A su vez, podemos escribir
f [x0 , x1 ] − f [x, x0 ]
f [x, x0, x1 ] =
x1 − x
que podemos reescribir como

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

y nos permite expresar f (x) como

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

Repitiendo este proceso nveces, llegamos a la expresión

f (x) = f [x0 ]+(x−x0 ) f [x0 , x1 ]+(x−x0 ) · · · (x−xn−1 ) f [x0 , x1 , . . . , xn ]+(x−x0 ) · · · (x−xn ) f [x0 , x1 , . . . , xn , x]

Los n primeros términos constituyen el polinomio interpolador, dado por la fórmula de Newton,
y por lo tanto, el último término es el término de error. Tenemos por consiguiente la relación:

f (n+1) (xm )
(x − x0 ) · · · (x − xn ) f [x0 , x1 , . . . , xn , x] = (x − x0 ) · · · (x − xn )
(n + 1)!

de donde deducimos la siguiente relación entre las diferencias divididas y las derivadas:

f (n+1) (xm )
f [x0 , x1 , . . . , xn , x] =
(n + 1)!

En el caso de puntos igualmente espaciados la anterior igualdad toma la forma:

&n+1 f0 f (n+1) (xm )

=
(n + 1)!hn+1 (n + 1)!

Luego podemos escribir la siguiente relación entre las derivadas y las diferencias finitas:

&n f0 = hn f (n) (xm )

donde xm es un punto del intervalo [x0 , xn ]. Cuando el intervalo [x0 , xn ] se hace infinitamente
pequeño, tenemos el siguiente límite cuando xn −→ x0 :

lı́m &n f0 −→ hn f (n) (x0 )

xn −→x0

que nos dice que las diferencias y las derivadas de orden n son proporcionales en este límite.
CAPÍTULO 5. INTERPOLACIÓN 65

5.5. Elección de los puntos de interpolación

El polinomio interpolador puede presentar comportamientos patológicos, en el sentido de
que, cuando aumentamos el número de puntos de interpolación en un intervalo, la desviación
máxima entre el polinomio y la función en dicho intervalo aumente con el grado del polinomio.
Esto pasa con las funciones de Bernstein, B(x) = |x|, y Runge

1
R(x) =
1 + 25x2
en el intervalo [−1, 1]. Cuando aumentamos el grado del polinomio interpolador de estas funcio-
nes con puntos igualmente espaciados, sucede que

máxn→∞ | f (x) − Pn (x)| −→ ∞

−1 ≤ x ≤ 1

Podemos preguntarnos cual es la mejor elección de los puntos de interpolación, en caso de

que tengamos libertad para elegirlos. La respuesta es que el error se minimiza si se eligen los
puntos de interpolación como los ceros de los polinomios de Chebychev, definidos por

Tn (x) = cos(n arc cos(x))

De su definición, es obvio que lox Tn (x) sólo toman valores en el intervalo [−1, 1], pues fuera
de este intervalo la función arc cos(x) no está definida. Como sus ceros están en el intervalo
[−1, 1], para una función f (x) definida en un intervalo [a, b] tenemos que realizar un cambio
b−a b+a
de variables a la variable t definida por x = t+ , de forma que f (t) toma valores en
2 2
[-1,1]. Los valores de los ceros de los polinomios de Chebychev vienen dados por los valore de
x que satisfacen
cos(nθ) = 0
donde θ = arc cos x. Tenemos por lo tanto

(2k + 1)π
nθ =
2
(2k + 1)π
θ=
2n
por lo que los ceros xn de Tn (x) vienen dados por la expresión:
& '
(2k + 1)π
xn = cos
2n

Los polinomios de Chebychev satisfacen la relación de recurrencia

Tn+1 (x) = 2xTn (x) − Tn−1 (x)

CAPÍTULO 5. INTERPOLACIÓN 66

con T0 (x) = 1 y T1 (x) = x. Esta relación se puede demostrar fácilmente a partir de las relaciones
trigonométricas

cos(n + 1)θ = cos nθ cos θ − sin nθ sin θ

cos(n − 1)θ = cos nθ cos θ + sin nθ sin θ

Sumando ambas igualdades, y poniendo θ = arc cos x, tenemos

cos((n + 1) arc cos x) = 2x cos(n arc cos x) − cos((n − 1) arc cos x)

que da la anterior relación de recurrencia, teniendo en cuenta la definicion de Tn (x). De la relación

de recurrencia se obtiene que el coeficiente del término principal de Tn (x) es 2n−1 .
La razón de la idoneidad de los ceros de los polinomios de Chebychev como puntos de
interpolación radica en que los polinomios de Chebychev son los polinomios con coeficiente
principal 2n−1 con un menor valor máximo en [-1,1]. Esta afirmación es fácil de demostrar. Los
polinomios de Chebychev alcanzan su valor máximo (+1 o -1), en x = 1, x = −1, y en los puntos
que anulan sus derivadas (T % (x) = 0). Estos puntos vienen dados por
kπ
xk = cos
n
con k = 1, 2,. . . , n−1, por lo que tenemos que alcanzan su máximo valor absoluto en n+1 puntos,
en los que alternan el signo. Si suponemos que un polinomio Qn (x), con el mismo coeficiente del
término principal, es siempre menor en módulo que Tn (x), entonces la diferencia Tn (x) − Qn (x),
que es un polinomio de grado n − 1 (los términos de grado n de ambos polinomios son iguales),
debe de tomar valores con signos opuestos en todos los máximos (en valor absoluto) de Tn (x), o
sea en n + 1 puntos. Como esta diferencia es una función continua, debe de anularse en n puntos.
Pero un polinomio de grado n − 1 que se anula en n puntos debe ser idénticamente nulo. Por lo
tanto, no existe un polinomio de coeficiente principal 2n−1 que sea menor que Tn (x) en [-1,1].
Consideremos ahora el término del error del polinomio interpolador
n
f (n+1) (α)
f (x) = Pn (x) + ∏ (x − x j )
(n + 1)!
j=0

El término ∏n (x − x j ) es Tn+1 (x)/2n cuando los x j se toman como los ceros de Tn (x), ya
j=0
que tiene los mismos ceros que Tn+1 (x) y debe coincidir con él, salvo una constante multiplicati-
va, y esta constante debe ser 2−n para que el coeficiente del término principal sea 1. Por lo tanto,
si se toman los puntos de interpolación en los ceros de los polinomios de Chebychev, el error de
Pn (x) viene dado por
Tn+1 (x) f (n+1) (α)
f (x) = Pn (x) +
2n (n + 1)!
Por lo tanto, el término de error tiene la propiedade de Tn+1 (x) de tomar el menor valor absoluto
posible, en [-1,1], si despreciamos las variaciones debidas a f (n+1) (α) (αdepende del punto x en
el que se calcula el error).
CAPÍTULO 5. INTERPOLACIÓN 67

Cuando se interpolan las funciones de Bernstein y Runge en los ceros de los polinomios de
Chebychev se obtiene que
máx | f (x) − Pn (x)| −→ 0
−1 ≤ x ≤ 1 n→∞
aunque la convergencia es lenta, en comparación con la que se obtiene para funciones de com-
portamiento “normal”.

4.6. Interpolación por splines

La interpolación polinomial tiene el inconveniente de que el polinomio interpolador puede
oscilar fuertemente entre los puntos interpolados, como hemos visto en el apartado anterior. Para
muchas aplicaciones prácticas interesa un algoritmo que proporcione una función interpoladora
que se comporte suavemente. El modelo más sencillo es una interpolación lineal a tramos, cons-
tituida por rectas que unen los puntos interpolados. Sin embargo, este tipo de aproximación tiene
una derivada discontinua en lo puntos de interpolación. Para muchas aplicaciones prácticas, en
las que se incluyen la resolución de ecuaciones diferenciales en las que intervienen funciones
medidas experimentalmente, interesa una función aproximadora que tenga derivadas continuas
hasta un orden dado, aparte de ser continua. El método de splines se basa en encontrar funciones
que cumplan estas características. Se define un spline de orden m en una serie de n + 1 puntos
de interpolación {(x0 , y0 ), . . . (xn , yn )} como un conjunto de n funciones Sk (x) (frecuentemente
polinomios de orden m) definidas en el intervalo [xk , xk+1 ] , que satisfacen

Sk (xk ) = Sk+1 (xk )

Sk% (xk ) = %
Sk+1 (xk )
Sk%% (xk ) = %%
Sk+1 (xk )
.. .. ..
. . .
(m−1) (m−1)
Sk (xk ) = Sk+1 (xk )

Como veremos más abajo, hay que introducir en general condiciones de contorno adicionales en
los extremos x0 y xn para que las funciones Sk (x) estén unívocamente definidas. Los splines más
utilizados son los cúbicos, debido a que son los más sencillos con derivada segunda continua, y
por lo tanto son adecuados para aproximar funciones que intervienen en ecuaciones diferenciales
de segundo orden. Si escribimos

Sk (x) = ak + bk (x − xk ) + ck (x − xk )2 + dk (x − xk )3

las dos primeras derivadas vienen dadas por

Sk% (x) = bk + 2ck (x − xk ) + 3dk (x − xk )2

Sk%% (x) = 2ck + 6dk (x − xk )

CAPÍTULO 5. INTERPOLACIÓN 68

La condición de que los splines pasen por los puntos de interpolación da Sk (xk ) = ak = yk , que nos
dice que los coeficientes ak son los valores de la función interpolada en los nodos. El significado
de los otros coeficientes es obviamente Sk% (xk ) = bk , S%% (xk ) = 2ck y S%%% (xk ) = 6dk . Tenemos que
las ecuaciones de continuidad de las funciones Sk (x) y sus derivadas primera y segunda en los
n − 1 puntos de interpolación intermedios xk+1 , k = 0, . . . , n − 2 (las condiciones de continuidad
no se aplican a los extremos) dan

ak+1 = ak + bk (xk+1 − xk ) + ck (xk+1 − xk )2 + dk (xk+1 − xk )3

bk+1 = bk + 2ck (xk+1 − xk ) + 3dk (xk+1 − xk )2
2ck+1 = 2ck + 6dk (xk+1 − xk )

que constituyen un sistema de 3(n − 1) ecuaciones para determinar los 3n coeficientes bk ,ck y dk .
Una condición adicional es que el último spline pase por el último punto: Sn−1 (xn ) = yn . Nece-
sitaremos por lo tanto dos condiciones adicionales para poder determinar todos los coeficientes.
Daremos más adelante estas condiciones. Para simplificar la notación, introducimos la definición
hk = xk+1 − xk . Podemos escribir las ecuaciones de continuidad anteriores como
ak+1 − ak
= bk + ck hk + dk h2k
hk
bk+1 − bk = 2ck hk + 3dk h2k
ck+1 − ck = 3dk hk (5.1)

De las posibles formas de resolver el sistema de ecuaciones anterior, la más comoda es despejar
los coeficientes dk y bk de forma que obtengamos un sistema de ecuaciones para los coeficientes
ck . Si despejamos
ck+1 − ck
dk =
3hk
en la última ecuación y lo introducimos en la segunda ecuación, obtenemos
ck+1 − ck 2
bk+1 − bk = 2ck hk + 3 hk = (ck + ck+1 )hk (5.2)
3hk
Necesitamos otra expresión independiente de bk+1 − bk para obtener una ecuación para los ck .
En la primera de las ecuaciones 5.1 podemos despejar bk :
ak+1 − ak ak+1 − ak ck+1 − ck 2 ak+1 − ak 2 1
bk = − ck hk − dk h2k = − ck hk − hk = − ck hk − ck+1 hk
hk hk 3hk hk 3 3
y cambiando k por k + 1 obtenemos
ak+2 − ak+1 2 1
bk+1 = − ck+1 hk+1 − ck+2 hk+1
hk+1 3 3
Restando estas dos ecuaciones queda
ak+2 − ak+1 ak+1 − ak 1 1 2 2
bk+1 − bk = − − ck+2 hk+1 + ck+1 hk − ck+1 hk+1 + ck hk (5.3)
hk+1 hk 3 3 3 3
CAPÍTULO 5. INTERPOLACIÓN 69

Igualanddo esta expresión con la ecuación 5.2, obtenemos finalmente una relación entre los ck :
1 2 1 ak+2 − ak+1 ak+1 − ak
ck+2 hk+1 + ck+1 (hk + hk+1 ) + ck hk = −
3 3 3 hk+1 hk
Estas ecuaciones forman un sistema tridiagonal. Si definimos
& '
ak+2 − ak+1 ak+1 − ak
vk+1 = 3 −
hk+1 hk
tenemos que el sistema tridiagonal queda como
ck+2 hk+1 + 2ck+1 (hk + hk+1 ) + ck hk = vk+1
cuya forma matricial es:
    
h0 2(h0 + h1 ) h1 0··· 0 0 0 c0 v1
 0 h1 2(h1 + h2 ) h2 0 ··· 0  c1   v2
    
 0 0 h2 2(h2 + h3 ) h3 ··· 0  c2   v3
    
 0 0 0 h3 2(h3 + h4 ) ··· 0  .. =  ..
  .    .
 .. .. .. .. .. .. ..    
 . . . . . . .   cn−2   vn−3 
0 0 0 ··· hn−3 2(hn−3 + hn−2 ) hn−2 cn−1 vn−2
(5.4)
que es un sistema de n − 2 ecuaciones para n incógnitas. Necesitamos dos ecuaciones adicionales
para determinar los ck . Estas dos ecuaciones se proporcionan mediante dos condiciones en cierta
forma arbitrarias, pero de las que depende la calidad del ajuste. Las dos formas más usuales de
proporcionar estas condiciones adicionales son:
1. La condición denominada splines naturales, que impone derivada nula en los extremos:
S0%% (x0 ) = Sn−1
%% (x ) = 0.
n

2. La condición de splines sujetos (clamped), que fija la derivada en los extremos S0% (x0 ) =
f % (x0 ) y Sn−1
% (xn ) = f % (xn ).
La condición de splines naturales implica c0 = cn = 0. El coeficiente cn no pertenece a ninguno
de los splines en [x0 , xn ], sino a un spline arbitrario que comienza en el extremo xn , pero la
introducción de Sn (x) permite extender las ecuaciones para los coeficientes ck con una ecuación
adicional
hn−2 cn−2 + 2(hn−2 + hn−1 )cn−1 = vn−1
La condición c0 = 0 nos permite eliminar la primera columna de la matriz del sistema 5.4. Con
estas dos modificaciones, el sistema de ecuaciones queda de la forma:
 
2(h0 + h1 ) h1 0··· 0 0 0    
  c1 v1
 h1 2(h1 + h2 ) h2 0 ··· 0    
  c2   v2 
 0 h2 2(h2 + h3 ) h3 ··· 0    
  c3   v3 
 0 0 h3 2(h3 + h4 ) ··· 0  .. = .. 
 .. .. .. .. .. ..  .   . 
 . . . . . .    
  cn−2   vn−2 
 0 0 ··· hn−3 2(hn−3 + hn−2 ) hn−2 
cn−1 vn−1
0 0 0 ··· hn−2 2(hn−2 + hn−1 )
CAPÍTULO 5. INTERPOLACIÓN 70

que es un sistema de n − 1 ecuaciones con n − 1 incógnitas. Este sistema se puede resolver

fácilmente por eliminación de Gauss. Cada paso de eliminación se aplica sólamente a la fila
inmediatamente inferior, ya que los elementos inferiores de la columna son nulos. El factor de
eliminación de la segunda fila es
h1
f2 = −
2(h0 + h1 )
Después de la eliminación, la subdiagonal inferior se anula. El término en la diagonal toma un
nuevo valor u2 dado por
h21
u2 = 2(h1 + h2 ) −
2(h0 + h1 )
Los elementos uk de la diagonal a partir de la segunda fila, después de la reducción, vienen
dados por
h2
uk+1 = 2(hk + hk+1 ) + hk fk+1 = 2(hk + hk+1 ) − k
uk
hk
fk+1 = −
uk
y los términos independientes quedan como,
hk vk
v%k+1 = vk+1 −
uk
El sistema de ecuaciones, después de la reducción de Gauss, toma la forma:
  
2(h0 + h1 ) h1 0 · · · 0 0 0 c1  
 v1
 0 u2 h2 0 ··· 0   
  c2   v% 
 0     2 
 0 0 u3 h3 · · ·   c.3   v% 
 0     3 
 0 0 0 u4 · · ·   ..  =  .. 
 .. .. .. .. .. ..     
 . . . . . .   ...   % . 
    vn−2 
 0 0 ··· 0 un−2 hn−2   cn−2 
v%n−1
0 0 0 ··· 0 un−1 cn−1

Este sistema se resuelve fácilmente por sustitución hacia atrás,

v%n−1
cn−1 =
un−1
v% − hk ck+1
ck = k
uk
con lo que determinamos todos los coeficientes ck . A partir de los ck , determinamos los coefi-
cientes dk mediante la relación
ck+1 − ck
dk =
3hk
CAPÍTULO 5. INTERPOLACIÓN 71

Los coeficientes bk se obtienen de la relación de recurrencia

bk+1 = bk + (ck + ck+1 )hk

con el valor inicial de b0 determinado por la ecuación de continuidad de S0 y S1 en el primer nodo

x1 :
a1 − a0 a1 − a0 c1
b0 = − d0 h20 = − h0
h0 h0 3
donde hemos tenido en cuenta que c0 = 0.
En el caso de splines sujetos, damos como datos adicionales los valores de f0% y fn% ( o estima-
ciones de los mismos). Tenemos por lo tanto S0% (x0 ) = b0 = f0% y Sn−1 (xn ) = bn−1 + 2cn−1 hn−1 +
3dn−1 h2n−1 = fn% . La ecuación de continuidad en x1 da

a1 − a0
c0 h0 + d0 h20 = − f0%
h0
La anterior expresión, introduciendo d0 en función de c0 y c1 , queda como
c1 − c0 2 a1 − a0
c0 h0 + h0 = − f0%
3h0 h0
& '
a1 − a0 %
(2c0 + c1 )h0 = 3 − f0
h0
Análogamente, en el punto xn imponemos la condición de que la derivada de Sn−1 (x) valga
f % (xn ):
bn−1 + 2cn−1 hn−1 + 3dn−1 h2n−1 = fn%
De la condición de que el spline Sn−1 pase por el extremo xn obtenemos
an − an−1
bn−1 + cn−1 hn−1 + dn−1 h2n−1 =
hn−1
Restando ambas ecuaciones eliminamos bn−1 con lo que obtenemos la relación,
an − an−1
cn−1 hn−1 + 2dn−1 h2n−1 = fn% −
hn−1
Poniendo dn−1 en función de los ck ,
cn − cn−1
dn−1 =
3hn−1
obtenemos finalmente la ecuación adicional para cn buscada:
& '
% an − an−1
(cn−1 + 2cn )hn−1 = 3 fn −
hn−1
CAPÍTULO 5. INTERPOLACIÓN 72

Añadiendo al sistema de ecuaciones que determina los ck las dos ecuaciones adicionales en los
extremos inicial y final, tenemos que los ck vienen dados por la solución de un sistema de n + 1
ecuaciones con n + 1 incógnitas:
    
2h0 h0 0··· 0 0 0 c0 v0
 h0 2(h0 + h1 ) h1 0 ··· 0     
   c1   v1 
 0 h1 2(h1 + h2 ) h2 ··· 0     
   c2   v2 
 .. .. .. .. .. ..   ..  =  .. 
 . . . . . .     
  .   . 
 0 0 ··· hn−2 2(hn−2 + hn−1 ) hn−1   cn−1   vn−1 
0 0 ··· 0 hn−1 2hn−1 cn vn
con & '
a1 − a0 %
v0 = 3 − f0
h0
& '
% an − an−1
vn = 3 fn −
hn−1
Para los otros valores de k, vk viene dado por la misma expresión que en el caso de splines
naturales:
& '
ak+2 − ak+1 ak+1 − ak
vk = 3 −
hk+1 hk
Por lo tanto, en el caso de splines sujetos obtenemos un sistema de n + 1 ecuaciones con n + 1
incógnitas, que se resuelve como en el caso de splines lnaturales. La primera etapa de eliminación
de Gauss da el sistema reducido
    
2h0 h0 0 · · · 0 0 0 c0 v0
 0 u1 h1 0 ··· 0  c1   v%1 
    
 0 0 u2 h2 · · · 0  c2   v%2 
    
 0 0 0 u3 · · · 0  .. = .. 
  .   . 
 .. .. .. .. .. ..    % 
 . . . . . .   cn−1   vn−1 
0 0 · · · 0 0 un cn v%n
con
3
u1 = 2h1 − h0
2
h 2
uk = 2(hk−1 + hk ) − k−1 k = 1, . . . , n − 1
uk−1
h2n−1
un = 2hn−1 −
un−1
hk−1 vk−1
v%k = vk − , k>0
uk−1
CAPÍTULO 5. INTERPOLACIÓN 73

Este sistema se resuelve por sustitución hacia atrás:

%
v
cn = n
un

v%k − hk ck+1
ck =
uk
con v%0 = v0 y u0 = 2h0 . Los coeficientes bk y dk se determinan como en el caso de splines
naturales.