Interpolación polinómica en una variable

Unisolvencia

Supongamos que el espacio vectorial en el que deseamos resolver un problema de interpolación está generado por cinco funciones, (Subscript[h, i]), 1 ≤ i ≤ 5, definidas por

Clear["Global`*"]

(Subscript[h, 1])[x_] := x - 5(Subscript[h, 2])[x_] := (x^3) - 5 x(Subscript[h, 3])[x_] := (x^5) - 5 (x^3)(Subscript[h, 4])[x_] := (x^7) - 5 (x^5)(Subscript[h, 5])[x_] := (x^9) - 5 (x^7)

Y los datos corresponden a los valores de la función en los puntos 0, 1, 2, 3 y 4. Entonces el determinante de Gram es

\(p = {0, 1, 2, 3, 4}; 
Det[Table[h\_i[p], {i, 5}]]\)

diferente de cero, con lo cual el problema es unisolvente. Sin embargo, si, en lugar de considerar como último dato el valor en el punto 4, consideramos el valor en un punto a, existen varios valores de a para los cuales el problema no es unisolvente. De antemano los valores que coinciden con los demás nodos. Pero hay más, veámoslo.

\(p = {0, 1, 2, 3, a}; 
polinomio = Det[Table[h\_i[p], {i, 5}]]\)

Roots[polinomio == 0, a]

Para que el problema de interpolación de datos lagrangianos sea unisolvente no puede usarse como quinto nodo ninguno de los cuatro primeros ni sus opuestos. Tampoco se pueden usar otros dos puntos más: ni (Sqrt[5]) = 2.236067977498705... ni  su opuesto.

Volviendo a los nodos 0, 1, 2, 3 y 4, vamos a interpolar valores de la función f(x) = (((x + 1))^4) en dichos puntos.

\(f[x_] = \((x + 1)\)\^4; 
p = {0, 1, 2, 3, 4}; \)

\(A = Transpose[Table[h\_i[p], {i, 5}]]; 
b = f[p]; 

    interpolante = Table[h\_i[x], {i, 5}].LinearSolve[A, b]\)

Plot[{interpolante, f[x]}, {x, 0, 4}, 
      
PlotStyle 
        -> {Dashing[{0.01, 0.02}], RGBColor[0, 0, 1]}];

Problemas de interpolación polinómica más habituales

La orden InterpolatingPolynomial permite obtener el polinomio de interpolación cuyos datos son lagrangianos, Hermite, Taylor o Hermite generalizado. Hay que dar los datos entre llaves y a continuación, sepada por coma, la variable del polinomio interpolador. Si en un nodo se tienen como datos el valor de la función y el de 5 derivadas consecutivas, por ejemplo, los 6 datos van en una llave que sigue al nodo. Veamos un ejemplo; vamos a calcular el polinomio cuyos datos en el punto 1 son {5,1,0}, es decir, vale 5, su derivada primera vale 1 y su derivada segunda vale 0; los datos en el punto 3 son {1,3} y, finalmente, su valor en 5 es -2. El polinomio lo damos en la variable t. Observe la forma diferente de escribir los datos, entre llaves,  cuando sólo hay un dato, el valor de la función, y cuando hay más de un dato.

\(datos = {{1, {5, 1, 0}}, {3, {1, 3}}, {5, \(-2\)}}; 

    pol[t_] = Expand[InterpolatingPolynomial[datos, t]]\)

Vamos a emplear esta orden para resolver directamente varios problemas clásicos de interpolación polinómica.

Interpolación polinómica lagrangiana

Polinomio que pasa por los puntos (1,1), (2,0), (4,0), (5,1) y (6,1.5):

\(datos1 = {{1, 1}, {2, 0}, {4, 0}, {5, 1}, {6, 1.5}}; 

    pol1[t_] = Expand[InterpolatingPolynomial[datos1, t]]\)

\(g1 = 
      ListPlot[datos1, PlotStyle -> {PointSize[0.03]}, 
        DisplayFunction -> Identity]; 

    g2 = Plot[pol1[t], {t, 1, 6}, \ DisplayFunction -> Identity]; \)

Show[{g1, g2}, DisplayFunction -> $DisplayFunction];

Polinomio que interpola los valores de la función Exp[x]/(6  (x^2) + 5) en los puntos -1, 0 ,1, 2 y 3:

\(datos2 = 
      Table[{x, E\^x\/\(6\ x\^2 + 5\)}, {x, \(-1.\), 3}]; 

    pol2[t_] = Expand[InterpolatingPolynomial[datos2, t]]\)

\(g3 = 
      ListPlot[datos2, PlotStyle -> {PointSize[0.02]}, 
        DisplayFunction -> Identity]; 

    g4 = Plot[pol2[t], {t, \(-3\), 4}, DisplayFunction -> Identity, 
        PlotStyle -> {Dashing[{0.01, 0.02}]}]; 

    g5 = Plot[E\^x\/\(6\ x\^2 + 5\), {x, \(-3.\), 4}, 
        DisplayFunction -> Identity]; \)

Show[{g3, g4, g5}, DisplayFunction -> $DisplayFunction];

Interpolación polinómica de Hermite clásica

Polinomio que pasa por los puntos (1,1), (2,0), (4,0), (5,1) y (6,1.5) y su derivada primera vale cero en todos ellos:

\(datos3
       = {{1, {1, 0}}, {2, {0, 0}}, {4, {0, 0}}, 
        
{5, {1, 0}}, {6, {1.5, 0}}}; 

    pol3[t_] = Expand[InterpolatingPolynomial[datos3, t]]\)

\(g6 = Plot[pol3[t], {t, 1, 6.1}, DisplayFunction -> Identity]; 

    Show[{g6, g1}, DisplayFunction -> $DisplayFunction]; \)

Polinomio que interpola los valores de la función Exp[x]/(6  (x^2) + 5) y los de su derivada en los puntos -1, 0 ,1, 2 y 3:

\(datos4 = 
      Table[{x, {E\^x\/\(6\ x\^2 + 5\), 
            ∂\_{z, 1}\(E\^z\/\(6\ z\^2 + 5\)\) /. 
              z -> x}}, {x, \(-1.\), 3}]\)pol4[t_] = Expand[InterpolatingPolynomial[datos4, t]]

g7 = 
      Plot[pol4[t], {t, -2, 4}, 
        PlotStyle -> {Dashing[{0.025, 0.025}]}, 
        DisplayFunction -> Identity];

Show[{g3, g5, g7}, DisplayFunction -> $DisplayFunction];

Se observa una mayor coincidencia de este interpolante y la función interpolada en el intervalo [-1, 3] que del interpolante de Lagrange en los mismos puntos, visto más arriba. Sin embargo, fuera de dicho intervalo la diferencia es mayor en este caso.

Interpolación polinómica de Taylor

Finalmente vamos a calcular el interpolante de Taylor, de grado 5, para la misma función que venimos considerando, y centrado en el punto x = 1 (nodo central de los considerados anteriormente).

\(datos5
       = {{1, Table[
            ∂\_{z, i}\(E\^z\/\(6\ z\^2 + 5\)\) /. 
              z -> 1., {i, 0, 4}]}}\)pol5[t_] = Expand[InterpolatingPolynomial[datos5, t]]

g8 = 
      Plot[pol5[t], {t, -1, 3.}, 
        PlotStyle -> {Dashing[{0.025, 0.025}]}, 
        DisplayFunction -> Identity];

Show[{g5, g8}, DisplayFunction -> $DisplayFunction];

Mathematica dispone de una orden para hallar el desarrollo de Taylor de una función hasta el término que se desee. su sintaxis es Series[función,{variable, origen, orden}]. Dicho desarrollo se efectúa según potencias de variable-origen, e incluye un término de error, que se puede eliminar incluyendo la orden Normal. Por ejemplo, el desarrollo de Taylor de la función seno en el punto cero de orden tres es

Series[Sin[x], {x, 0, 3}]

El correspondiente polinomio de Taylor es

Normal[%]

Fórmulas de Lagrange y Newton para la interpolación polinómica

Vamos a comparar un aspecto importante entre la fórmula de Lagrange y la de Newton para la interpolación polinómica lagrangiana: el tiempo de evaluación del polinomio obtenido mediante ambas representaciones. Para ello basta comparar los tiempos que tarda el ordenador en dibujar sus gráficas correspondientes. Lógicamente los polinomios se dejan escritos tal como aparecen en las fórmulas, es decir, sin pasarlos a la base canónica.
Los nodos de interpolación van a ser 0, 1, 2, ..., 10, y los guardaremos en una lista x. Los valores a interpolar serán números aleatorios comprendidos entre 0 y 8, que guardamos en un vector y. Se generan a partir de la orden Random[].

Clear["Global`*"]

\(n = 12; 
x = Table[i, {i, 0, n}]; 

    y = Table[8\ Random[], {i, 0, n}]\)

El polinomio de interpolación de Lagrange se puede formar directamente:

pLagrange[t_] := 
      (∑_(i = 1)^(n + 1))y[[i]] 
          ( (∏_(j = 1)^(i - 1))((t - 
                  x[[j]])/(x
                    [[i]] - 
                  x[[j]]))) 
          ( (∏_(j = i + 1)^(n + 1))((t - 
                  x[[j]])/(x
                    [[i]] - 
                  x[[j]])));

Para comprobar que la interpolación ha sido correcta podemos evaluar la función precedente  en la lista de nodos, x, debiendo obtener el vector de valores interpolados, y.

pLagrange[x]

Podemos representar la gráfica del polinomio de interpolación que hemos calculado y medir el tiempo que se necesita para obtenerla, que dependerá, lógicamente, de las prestaciones del equipo que estemos empleando. Lo interesante será comparar, para un mismo equipo, los tiempos empleados cuando se utilizan las fórmulas de Lagrange y de Newton. La orden correspondiente es Timing, que produce una lista de salida con dos posiciones. Tenga paciencia.

Timing[
      Plot[pLagrange[t], {t, x[[1]], 
          x[[n + 1]]}]]

El núcleo de Mathematica ha empleado un tiempo elevado para producir la gráfica. Más adelante comentaremos la razón.

Vamos a resolver el mismo problema de interpolación usando ahora la fórmula de Newton. Para ello, hay que obtener primero los coeficientes de cada uno de los sumandos de la expresión de dicho polinomio, las llamadas diferencias divididas, que serán guardadas en la diagonal de una matriz, A, que a tal efecto crearemos. La primera columna de A será el vector y; las restantes serán copias de la primera, por ejemplo.

\(A = Table[y, {i, n + 1}]; 
A = Transpose[A]; \)

El siguiente procedimiento produce como resultado las diferencias divididas.

For[i = 2, i ≤ n + 1, i++, 
      
      
      For[j = i, j ≤ n + 1, j++, 
        
             
        A[[j, i]] = 
          ((A[[j, i - 1]] - 
              A[[j - 1, i - 1]])/(x
                [[j]] - 
              x[[j - i + 1]]))]]

Ahora podemos formar el polinomio de interpolación haciendo uso de la fórmula de Newton.

pNewton[t_] := 
      (∑_(i = 1)^(n + 1))A[[i, i
            ]] 
          (∏_(j = 1)^(i - 1))(
              t - x[[j]]);

Para ver si los dos métodos descritos conducen al mismo polinomio, simplificamos la
diferencia de los polinomios calculados.

Simplify[pNewton[t] - pLagrange[t]]

Los errores de redondeo hacen que no obtengamos el polinomio cero, pero los coeficientes involucrados tienen unas magnitudes muy pequeñas. De hecho, si evaluamos la expresión anterior con Chop se hacen cero los del orden de (10^(-16)).

Chop[%]

Dibujemos el polinomio de interpolación calculado a partir de la fórmula de Newton pidiendo información del tiempo necesario.

Timing[
      Plot[pNewton[t], {t, x[[1]], 
          x[[n + 1]]}]]

La representación gráfica ha requerido un tiempo muy inferior al que necesitamos cuando el polinomio estaba escrito con la fórmula de Lagrange. Ello se debe a que los sumandos de la representación de Newton son polinomios de grados inferiores (excepto el último) a los que intervienen en la fórmula de Lagrange, lo que redunda en un tiempo de cálculo inferior.

Comportamiento asintótico de la interpolación polinómica. Ejemplo de Runge

¿Qué ocurre cuando se aumenta el número de nodos de interpolación? 
Supongamos que en el intervalo [-5,5] elegimos nodos uniformemente distribuidos en un número creciente y que interpolamos la función f(x) = (1/((x^2) + 1)) en esos nodos. Vamos a medir la distancia entre esta función y su polinomio de interpolación en el punto medio de los dos primeros nodos. Haremos uso de la expresión de Newton dado su  menor coste computacional.

Comenzaremos indicando el procedimiento que debe seguirse para representar el polinomio de interpolación de grado menor o igual que 10, por ejemplo, y luego calculamos los valores en el punto medio de los dos primeros nodos de los errores de interpolación asociados a los polinomos de grados 10, 20, ... , 80 , a fin de apreciar el comportamiento.

Definimos la función.

f[z_] := (1/((z^2) + 1))

Formamos las listas de nodos (uniformemente espaciados en [-5,5]) y de valores de f(z) que deben ser interpolados. En ellos calculamos las diferencias divididas y formamos el polinomio de interpolación.

\(n = 10; 
h = 10\/n; 
x = Table[\(-5.\) + i\ h, {i, 0, n}]; 

    y = f[x]; 
A = Table[y, {i, n + 1}]; 
A = Transpose[A]; 

    For[i = 2, i ≤ n + 1, \(i++\), 
      
\ \ \ \ \ \ \ 
      For[j = i, j ≤ n + 1, \(j++\), 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        A[[j, i]] = 
          \(A[[j, i - 1]] - 
              A[[j - 1, i - 1]]\)\/\(x
                [[j]] - 
              x[[j - i + 1]]\)]]\)pNewton[t_] = 
      (∑_(i = 1)^(n + 1))A[[i, i
            ]] 
          (∏_(j = 1)^(i - 1))(
              t - x[[j]]);

Representamos la función interpolada y el interpolante.

Plot[{f[t], pNewton[t]}, {t, -5, 5}, 
      
PlotStyle 
        -> {{RGBColor[1, 0, 0], Thickness[0.008]}, 
          
                       {
            RGBColor[0, 0, 1]}}];

La aproximación es buena en un amplio intervalo alrededor del origen, como muestra la gráfica anterior. Sin embargo, el comportamiento del aproximante cerca de los extremos es muy deficiente (observe la siguiente gráfica, relativa al error: f-pNewton).

Plot[f[t] - pNewton[t], {t, -5, 5}, 
      PlotStyle -> {RGBColor[1, 0, 0]}];

Para hacernos una idea de la magnitud de los errores anunciados, vamos a calcular para los diferentes grados 10, 20, ..., 80 el error en el punto medio entre los dos primeros nodos.

For[n = 10, n ≤ 80, n = n + 10, 
      
       h = (10/n); 
      
      x = Table[-5. + i h, {i, 0, n}]; 
      
     y = f[x]; 
     
      A = Table[y, {i, n + 1}]; A = Transpose[A]; 
    
      For[i = 2, i ≤ n + 1, i++, 
        
           
        For[j = i, j ≤ n + 1, j++, 
          
                  
          A[[j, i]] = 
            ((A[[j, i - 1]] - 
                A[[j - 1, i - 1]])/(x
                  [[j]] - 
                x[[j - i + 1]]))]]; 
      
   
      p = (1/2) (
            x[[1]] + 
              x[[2]]); 
      
   
      error = f[p] - 
          (∑_(i = 1)^(n + 1))A[[i, i
                ]] 
              (∏_(j = 1)^(i - 1))(
                  p - x[[j]]); 
      
   
      Print["Con ", n + 1, " nodos el error en ", p, " es ", 
        error]]

Las salidas anteriores sugieren que el error crece sin límite, hecho probado por Runge para este ejemplo. Ello nos indica que no debemos utilizar la interpolación polinómica lagrangiana confiados en la convergencia con tal que la diferencia entre nodos tienda a cero.

Nodos de Chebyshev

En el error de la interpolación polinómica lagrangiana para una función suficientemente suave intervienen dos elementos: de un lado, la norma del máximo de una cierta derivada de la función, y la norma del polinomio obtenido al multiplicar los binomios t-x[[i]] y dividir por el factorial adecuado, del otro. Ésta última se minimiza cuando los nodos de interpolación son los ceros del correspondiente polinomio ortogonal de Chebyshev.

Los ceros de los polinomios de Chebyshev de grado mayor o igual que cero, tomados como nodos de interpolación, permiten construir polinomios de interpolación de grado creciente. La continuidad de la primera derivada de la función interpolada garantiza la convergencia a cero del error de interpolación. Sin embargo, la continuidad no es, en general, suficiente.

Generemos, por ejemplo, los 11 ceros del polinomio de Chebyshev de grado 11 correspondientes al intervalo [-5,5].

\(a = \(-5.\); b = 5.; 
m = 11; 

    xChebyshev = 
      N[Table[\(a + b\)\/2 + 
            1\/2\ \((a - b)\)\ Cos[\(\((2\ i - 1)\)\ π\)\/\(2\ m\)], {i, 
            m}]]\)

Representemos el polinomio de interpolación de f en estos puntos, para lo que vamos a hacer uso de la orden InterpolatingPolynomial.

datosChebychev = 
      Table[{xChebyshev[[i]], 
          f[xChebyshev[[i]]]}, {i, m}]pChebyshev = InterpolatingPolynomial[datosChebychev, t];

Plot[{f[t], pChebyshev}, {t, a, b}, 
      
        
      PlotStyle 
        -> {{RGBColor[1, 0, 0], Thickness[0.008]}, 
          
                                {RGBColor[0, 0, 1]}}];

Se observa que no se producen las grandes oscilaciones que surgían al considerar los nodos igualmente espaciados.

Podemos ilustrar la situación examinando los términos polinómicos del error correspondientes a los nodos de Chebyshev y a los uniformemente espaciados.

El primero se define mediante la siguiente orden.

wChebyshev[t_] = 
      (((∏_(i = 1)^m)(
            t - xChebyshev[[i]])
            )/(m!))

grafChebyshev = 
      Plot[wChebyshev[t], {t, a, b}, 
        
                                        PlotStyle -> {RGBColor[1, 0, 0]}];

Las oscilaciones se mantienen entre unos valores notablemente bajos y se comportan de manera regular.

Para el segundo,  formamos una lista con los nodos resultantes al subdividir el intervalo [-5,5] en 10 partes iguales (obtenemos 11 puntos), construimos el correspondiente producto de binomios y lo representamos.

\(n = 10; 

    xUniforme = Table[a + \(i\ \((b - a)\)\)\/n, {i, 0, n}]\/\(n!\); 

    wUniforme[t_] = 
      ∏\+\(i = 1\)\%\(n + 1\)\((
          t - xUniforme[[i]])\); 

    grafUniforme = 
      Plot[wUniforme[t], {t, a, b}, 
        
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {RGBColor[0, 0, 1]}]; \)

Es evidente el mal comportamiento cuando nos aproximamos a los extremos, incluso cuando estamos lejos de los mismos.

Show[grafChebyshev, grafUniforme];

Ejercicios

1.- Se considera la función no negativa Exp[-(3.9 + (5/2000)) t] (t^2). Calcule su polinomio de interpolación en los nodos 0,1,2,3,4,5. Represéntelo gráficamente y calcule su valor en el punto 4.75.

2.- Programe el cálculo de las diferencias divididas usando sólo el vector de valores de la función y otro de igual longitud.

3.- Interpole la función Exp[x] + (x/5) en los puntos -n,...,-1,0,1,...,n para n=1, 2, 3. Represente la función y los interpolantes en el intervalo [-3,3] utilizando diferentes colores.

4.- Para la misma función del ejercicio anterior, halle el interpolante de Hermite en el espacio de los polinomios de grado menor o igual que trece para los nodos -3,-2,-1,0,1,2 y 3. Calcule también el de Taylor en el mismo espacio con nodo x=0. Represéntelos en [-3,3]. Observe que la función de partida es creciente y los interpolantes no.

5.- Para qué valores del parámetro a es unisolvente el problema consistente en interpolar en el conjunto de los polinomios de grado 5 los datos g(0),g''(0), g'''(0), g'(1),g'(a) y g''(a) (Hermite-Birkhoff)

6.- Para qué valores del parámetro a es unisolvente el problema consistente en interpolar en el conjunto de los polinomios de grado 5 los datos g(0),g''(0), g'''''(0), g'(1),g'(a) y g''(a).