Uno de los principales usos de la interpolación ha sido el de hallar valores intermedios a los calculados en tablas trigonométricas
o astronómicas. Tal como dice por ejemplo el  Anuario del Observatorio Astronómico de 2003: "Muchas de las tablas de este
Anuario contienen listas de valores correspondientes a posiciones dadas para instantes de tiempo sucesivos de una duración de
un día. Por medio de la interpolación es posible determinar los valores de tales magnitudes para instantes intermedios a los que
aparecen en la tabla. El tipo de interpolación que hay que llevar a cabo depende de la precisión que se desea alcanzar, siendo quizás
el caso más difícil de tratar el del movimiento de la Luna".

    Para probar la existencia de este polinomio se pueden considerar sus coeficientes  (a_{0 ,}a_{1 ,}a_{2 ,}.... ,a_n ) como incógnitas
a determinar en las n+1 ecuaciones  p_n(x_i) = y_i , para n=0,1,2, ... ,n,   ya que los  (x_i,  y_i ) son números ya conocidos como
datos del problema. Aparece pues el sistema de  n+1 ecuaciones con  n+1 incógnitas
 

a0 + a1 x0  + a2 x02 + .... +an x0n  =  y0
a0 + a1 x1  + a2 x12 + .... +an x1n  =  y1
...
a0 + a1 xn  + a2 xn2 + .... +an xnn  =  yn

      Si llamamos  A   la matriz de este sistema, sabemos que el sistema tiene solución única si el determinante de A es distinto
de cero. Pero esta matriz es la conocida matriz de Vandermonde, cuyo determinante es distinto de cero si los x_i son distintos.

    Esta misma deducción nos puede servir como un método para el cálculo del polinomio de interpolación. Para hacerlo con
matlab hay un pequeño problema ya que no admite índices 0 en vectores y matrices. Debemos considerar pues un cambio de
notación en nuestros programas. Llamemos por ejemplo

al polinomio de interpolación que estamos buscando. En notación de matlab esto sería

                                        p(x) = c(1) * x^(n-1) + c(2) * x^(n-2) + ... + c(n-1) * x + c(n)

      Luego,  para hallar el vector de coeficientes c=(c₁,c₂, ..., c_n) ,  si  y=(y₁,y₂, ..., y_n) son los n datos. ( !Cuidado con el cambio
de notación!). Se tendrá
                                                           c' = A \ y'
esto es de A c^t = y ^t, c^t = A^-1 y ^t, donde A es la matriz de Vandermonde, donde ahora si  x = (x₁,x₂, ..., x_n)  son las abscisas de los
puntos a interpolar, será  A= ( a_ij ) =  (x_i^n-j). Para A(i,j)= a_ij , i representa fila, j columna.

    El programita siguiente calcula los coeficientes del polinomio de interpolación de y=sin x, tomando nueve nodos equidistantes
en el intervalo [-4,4]. El resultado aparece en la columna de la derechas, son los coeficientes del polinomio:
 

    Para evaluar el polinomio de interpolación p(x) en un punto t podemos utilizar lo que se llama método de Horner, que
consiste simplemente en poner
 

px = c(1);  for k = 2:n;   px = px .* t + c(k);  end

donde t  puede ser tanto un escalar como un vector de nodos. En matlab hay una función que hace lo mismo, px = polyval ( c, x ).
Luego si programamos en matlab podemos sustituir el código anterior por esta expresión. En el programa de abajo utilizamos
estos hechos para representar la gráfica del polinomio que interpola  9 nodos de y=sin(x) (los que están marcados con
asterisco).
 

n=9;  x=[-4:1:4];  y=sin(x);  plot(x,y,'*r')  hold on  for i = 1:n         for j = 1:n           A(i,j)=x(i)^(n-j);          end      end      c = A\y'      t=[-6:0.01:6];  px = c(1);  for k = 2:n;    px = px .* t + c(k);  end  plot(t,px)

   El polinomio de interpolación en forma de Newton

    Utilizar la matriz de Vandermonde para muchos nodos no es muy buena idea ya que el tiempo de cálculo para matrices grandes
es excesivo. Es mucho más sencillo utilizar el método clásico de las diferencias divididas de Newton. Recordemos su definición, para dos
nodos, se llama diferencia dividida de oreden uno a :

    Mientras que la diferencia dividida de orden n se obtiene por recurrencia a partir de las anteriores como

    El polinomio de Newton en differencias divididas es entonces
 

p(x)=f[x0]+(x-x0) f[x0,x1]+ (x-x0)(x-x1) f[x0,x1]+ +(x-x0)(x-x1) (x-xn-1) f[x0,x1, ... , xn]

    Las diferencias divididas son muy útiles para calcular el polinomio de interpolación. Para ello se suelen ordenasr en forma de tabla, tal como hacemos
aquí abajo, donce hemos mostrado el caso de tomar cuatro nodos. Los coeficientes del polinomio serían entonces los elementos de la diagonal superior.

x0   f[x0]  f[x0,x1]  x1  f[x1] f[x0,x1,x2]  f[x1,x2] f[x0,x1,x2,x3] x2  f[x2] f[x1,x2,x3] f[x2,x3] x3  f[x3] Tabla de diferencias divididas

    El programa siguiente calcula las diferencias divididas. El resultado puede verse a la derecha. Los coeficientes del polinomio
que hay que poner en la forma de Newton serían los elementos de la columna vertical  izquierda.
 

x=[-4:2:4]  y=sin(x)  for i = 1:n-1    for j=n:-1:i+1       y(j) = (y(j) - y(j-1)) / (x(j) - x(j-i));       fprintf('  %10.4f',y(j));    end    fprintf('\n');  end

x =       -4    -2     0     2     4  y =      0.7568   -0.9093         0    0.9093   -0.7568      -0.8330      0.4546      0.4546     -0.8330      -0.3219      0.0000      0.3219      -0.0537     -0.0537       0.0000

     Para practicar con estos métodos vamos a intentar buscar una curva famosa, una de las curvas que estudia Newton en
su clasificación de cúbicas, concretamente la llamada serpentina o anguinea
 

La serpentina de Newton es una cúbica de ecuación     Después de clasificar su 72ª curva, acaba Newton   diciendo:    "Y así son las especies en todo setenta y dos".  I. Newton Enumeratio Linearum Tertii Ordinis  (de la edición de De Analysi)

    Construimos el siguiente programa, que marca con asterisco unos puntos sobre la curva para luego hallar el polinomio de
interpolación. La curva de color azul es la gráfica de este polinomio, mientras que el otro trazo corresponde a la anguinea.
Utiliza la función newtoninterp(x,y,t), que se encuentra en el fichero newtoninterp.m
 

x=[-5:2:5]; % Nodos para interpolar  y=x./(1+x.^2);% Valores de la función en estos nodos  plot(x,y,'*')%Marca los puntos (xi,yi) a interpolar  hold on  t=[-6:0.01:6];% Puntos para dibujar  s=t./(1+t.^2);% Calcula la función exacta para dibujar  plot(t,s,'m')% Dibuja la función f(x)  axis([-6 6 -2.6 2.6]);  hold on  z=newtoninterp(x,y,t);% polinomio de interpolación  plot(t,z) %Gráfica del polinomio  legend('Puntos a interpolar','f = x/(1+x^2)','Polinomio de interp grado 5')  grid on

Problemas con la interpolación

    Hasta ahora hemos estudiado el polinomio de interpolación y parece que todo va bien. Una pregunta natural en este
contexto es la siguiente:

    Supongamos que dado un intervalo [a,b] lo vamos subdividiendo en más y más puntos, más concretamente tomemos

                         x_j=a+jh    para j=0,1,2,...,n    donde   h=(b-a)/n

y supongamos que construímos con estos puntos el polinomio de interpolación p_n(x) para una función dada f, esto es, que
p_n(x_i)=f(x_i),  para estos n puntos.
 
 

Supongamos que construímos el polinomio de interpolación pn(x) con estos puntos,      ¿Se tiene           cuando n tiende a infinito?

     La respuesta es negativa. C. Runge, en 1901, mostró un ejemplo donde esto no ocurre. Dada la función

        en el intervalo [-5,5]

si construímos el polinomio de interpolación p_n(x) en este intervalo, entonces seguro que no hay convergencia en los puntos
donde |x|>3.63. El problema parece pues que se tuerce, pero por otro lado se vuelve más interesante. Resulta que para ciertas
funciones, por ejemplo para f(x)=exp(x), si hay convergencia. Pero  para otras no. El problema con estas últimas es que
sus derivadas van creciendo demasiado en el intervalo considerado, esto es lo que sucede con esta función de Runge, que
parecía al principio bastante inofensiva. En las figuras que siguen mostramos este comportamiento:
 

Ejemplo de Runge. Aproximamos la función    por el polinomio de interpolación que pasa   por los puntos en asterisco. Se puede observar    que la aproximación es mala cerca de los extremos    del intervalo.

        Ya que la aproximación es mala en los extremos del intervalo, una idea para mejorar la aproximación es la
de olvidarnos de tomar nodos igualmente espaciados y tomar nodos que se concentren más cerca de los
extremos. De este modo al obligar al polinomio p_n(x) a pasar por estos puntos quizás se mejore la
aproximación. Por supuesto que tiene que haber un equilibrio en la disposición de los nodos  x_{i ,} pues si ponemos
pocos puntos en la región central del intervalo quizás perderíamos allí. Estas ideas son las que llevan
a una teoría de aproximación muy bonita, donde resulta que  los nodos a usar son los ceros de los
llamados polinomios de Chebyshev T_n(x).
 

La figura corresponde a la aproximación de la   misma  función de Runge anterior donde ahora   tomamos los nodos de  Chebyshev y la aproximación   se hace en el intervalo [-1 1] (para mayor comodidad   ya que los ceros de Chebyshev se definen en este  intervalo y no hay que hacer cambio de variable).    Obsérvese que los nodos están más concentrados   hacia los extremos del intervalo.

        P. L. Chebyshev  (1821-1894) ha sido uno de los mejores matemáticos rusos, conocido en gran parte por sus trabajos en
probabilidad, teoría de números, integración en un número de términos finitos y teoría de la aproximación. Fue gran amigo de
Liouville y de Hermite. Su interés en la teoría de la aproximación nació porque Chebyshev que era profesor de matemáticas en
San Petesburgo, fue profesor allí durante 35 años, viajó en 1852 a Paris y a Inglaterra donde estuvo muy interesado en la máquina de
vapor. Poco antes, James Watt había mejorado mucho el mecanismo de los pistones de la máquina. El mecanismo funciona
principalemente porque el pistón, que tiene un movimiento de vaivén rectilíneo, se tiene que transformar luego en circular,
porque la máquina al final tiene que mover una rueda. Esto se consigue por una serie de varillas articuladas (linkages) que
convierten movimiento rectilíneo en circular. Sin embargo, a pesar de que Watt había conseguido una buena aproximación,
se presentaban numerosos problemas. Esto era principalmente porque si miramos el movimiento al revés, es decir, de circular
al girar una rueda que se transforme mediante una serie de brazos articulados en rectilíneo, el movimiento que resulta
nos es completamente lineal, sino que presenta oscilaciones (curva de Watt), ello hacía que se resintiera la máquina.

      Chebyshev se dedicó de lleno a inventar un mecanismo para este propósito, que transformara movimiento rectilíneo
en circular. Encontró varios mecanismos aproximados, pero ninguno exacto. Varios años más tarde, en 1870, un joven estudiante
de Chebyshev llamado Lipkin, encontró un sencillo mecanismo que consigue el movimiento. Este mecanismo lleva también el
nombre del ingeniero francés Peaucellier, que encontró el mismo mecanismo en 1873. Es curioso observar que Chebyshev
no consiguiera encontrar la solución que buscaba con tanto esfuerzo, pero desarrolló en el camino unas preciosas matemáticas
que han tenido aplicación en otros problemas. Como el de interpolación que nos ocupa.
 
 

Mecanismo de Peaucellier, está formado por una serie de  brazos articulados, donde C está fijo y los demás son   pivotes movibles. Cuando el punto F' describe una línea,  el otro punto F  traza una circunferencia.

 

El siguiente ejemplo utiliza la función ginput de matlab, que toma puntos de la pantalla con el ratón. Por ejemplo, si usamos el código  xy=[ ];  x=[0:0.01:2];  plot(x,sin(x))  hold on  but = 1;      while but == 1          [xi,yi,but] = ginput(1);          plot(xi,yi,'ro')          hold on          n = n+1;          xy(:,n) = [xi;yi];      end dibuja y=sin(x), para luego ir tomando puntos sobre ella que va  mostrando en la consola.      La figura de la izquierda se ha hecho con una pequeña modificación de runge.m, que a su vez necesita las funciones difdivididas.m  y horner.m . El programa marca puntos sobre  la pantalla para a la vez trazar el polinomio de interpolación.

axis([0 10 0 10])      hold on      % Initially, the list of points is empty.      xy = [];      n = 0;      % Loop, picking up the points.      disp('Left mouse button picks points.')      disp('Right mouse button picks last point.')      but = 1;      while but == 1          [xi,yi,but] = ginput(1);          plot(xi,yi,'ro')          n = n+1;          xy(:,n) = [xi;yi];      end      % Interpolate with a spline curve.      t = 1:n;      ts = 1: 0.1: n;      xys = spline(t,xy,ts);      % Plot the interpolated curve.      plot(xys(1,:),xys(2,:),'b-');      hold off

Otro tipo de aproximación polinómica muy   usada en gráficos por ordenador es la spline    cúbica. Consiste en trazar una serie de puntos    que uniremos por pedazos de cúbica. Esto es,    tomamos una cúbica distinta que une cada par    de puntos consecutivos a intepolar. Los    coeficientes de cada cúbica se   tienen que tomar adecuadamente para que    hasta las segundas derivadas coincidan en    los puntos de enganche. El resultado es una   curva suave agradable a la vista.