Interpolación

    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".



    Problema de interpolación:
 
       Dados n+1 puntos  (x0 ,y0 ), (x1 ,y1 ), (x2 ,y2 ), ...., (xn ,yn ) del plano, hallar un polinomio de grado n

                            pn(x)= a0 + a1 x  + a2 x2 + .... +an xn

        que pase por estos puntos, esto es que pn(xi) = yi , para n=0,1,2, ... ,n

    Para probar la existencia de este polinomio se pueden considerar sus coeficientes  (a0 ,a1 ,a2 ,.... ,an ) como incógnitas
a determinar en las n+1 ecuaciones  pn(xi) = yi , para n=0,1,2, ... ,n,   ya que los  (xi,  yi ) 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 x0 =  y0
           a0 + a1 x1  + a2 x12 + .... +an x1 =  y1
           ...
           a0 + a1 xn  + a2 xn2 + .... +an xn =  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 xi 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=(c1,c2, ..., cn) ,  si  y=(y1,y2, ..., yn) son los n datos. ( !Cuidado con el cambio
de notación!). Se tendrá
                                                           c' = A \ y'
esto es de A ct = y t, ct = A-1 y t, donde A es la matriz de Vandermonde, donde ahora si  x = (x1,x2, ..., xn)  son las abscisas de los
puntos a interpolar, será  A= ( aij ) =  (xin-j). Para A(i,j)= aij , 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:
 

 n=9;
 x=[-4:1:4];
 y=sin(x);
 plot(x,y,'*')
 hold on
 for i = 1:n
        for j = 1:n
          A(i,j)=x(i)^(n-j); 
        end
     end
     c = A \ y'
c =

   -0.0000 
   -0.0001
    0.0000
    0.0077
   -0.0000
   -0.1649
    0.0000
    0.9988
         0

    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.
x  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

                         xj=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 pn(x) para una función dada f, esto es, que
pn(xi)=f(xi),  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 pn(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 pn(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  xi , 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 Tn(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.



Un magnífico programa: ginterp.m (de D. Arnold). Una aplicación java: http://www.wam.umd.edu/~petersd/interp.html

         Tomeu Barceló -- Departamento de Matemáticas -- Universidad Autónoma de Madrid