Interpolación spline

Introducción de los datos del problema

Especificamos unos nodos concretos en orden creciente (n en total, que darán lugar a n-1 subintervalos de longitudes (Subscript[h, k]) ,  k = 1, 2, ..., n - 1)

\(x = {0, π\/4, π\/2, \(3  π\)\/4, π} // N
\)\(n = Length[x]
\)\(a = x[\([1]\)]; \ b = x[\([n]\)]; 


    h = Table[x[\([k + 1]\)] - x[\([k]\)], {k, 1, n - 1}]\)

O bien consideramos una partición uniforme de cierto intervalo [a,b], con n-1 subintervalos y n nodos en total, incluidos los extremos del intervalo. Ahora las longitudes de los subintervalos serán todas iguales de longitudes (Subscript[h, k])  = ((b - a)/(n - 1)),  k = 1, 2, ..., 
      n - 1

a = 0.;  b = 2  π;

\(n = 20; 
x = Table[a + i\ \(b - a\)\/\(n - 1\), {i, 0, n - 1}]\)h = Table[x[[k + 1]] - x[[k]], {k, 1, n - 1}]

Y unos datos, que pueden ser los valores que una función concreta toma en los nodos indicados.

\(Clear[f, z]
\)\(f[z_] := Sin[z]
\)y = f[x] // Chop

Y cuya gráfica podemos guardar para posteriores comparativas

graf = 
      Plot[f[t], {t, x[[1]], 
          x[[n]]}, 
        PlotStyle -> {RGBColor[1, 0, 0]}];

También podríamos guardar tan sólo los valores puntuales que hemos considerado, en ese caso, o bien introducirlos nosotros mismos en la lista y={(Subscript[y, 1]), ...,  (Subscript[y, n - 1]),  (Subscript[y, n])}.

grafpuntos = 
      ListPlot[Table[{x[[k]], y[[k]]}, {k, n}], 
        PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.02]}];

Utilizamos potencias truncadas para hallar el spline de interpolación. Se definen, por ejemplo, del siguiente modo:

potenctruncada[z_, m_] := Which[z < 0, 0, True, (z^m)]

Podemos trasladar al punto x0∈ℝ arbitrario muy fácilmente:

Plot[potenctruncada[z - 2, 3], {z, -5, 5}];

Plot[potenctruncada[z - 10, 3], {z, 9, 12}, 
      AspectRatio -> Automatic];

Interpolación spline de Lagrange

Spline lineal de interpolación (clase (C^0))

En esta sección se muestra cómo obtener un spline lineal de interpolación una vez especificados los nodos de interpolación, (Subscript[x, 1]) < (Subscript[x, 2.]).. < (Subscript[x, n - 1]) < (Subscript[x, n]), y los valores que deben ser interpolados,  (Subscript[y, 1]), (Subscript[y, 2.])..,  (Subscript[y, n - 1]), (Subscript[y, n]).

Como la restricción del spline de interpolación a cada uno de los subintervalos [(Subscript[x, k]), (Subscript[x, k + 1])] es un polinomio de grado menor o igual que uno, la determinamos directamente a partir de los valores en los extremos.

En consecuencia, debemos considerar un total de n-1 subintervalos.

Usando potencias truncadas es muy fácil definir el spline. Comenzamos escribiéndolo como combinación lineal de los elementos de la base apropiada:

spline = 
      (∑_(i = 0)^1) c[i] (z^i) + 
        (∑_(i = 2)^(n - 1))c[i] potenctruncada[z - x[[i]], 1];

Se imponen las condiciones de interpolación, obteniendo un sistemas de ecuaciones en las incógnitas consideradas:

sistema = 
      Table[(spline /. z -> x[[i]]) == y[[i]], {i, n}]

Las incógnitas son

incognitas = Table[c[i], {i, 0, n - 1}]

Resolvemos el sistema:

sol=Solve[sistema,incognitas]

Representamos gráficamente la solución del problema de interpolación.

grafspl1 = Plot[spline /. sol[[1]], {z, a, b}];

Y la mostramos conjuntamente con la función interpolada.

Show[graf, grafspl1];

Podemos representar muy fácilmente el error cometido.

graferror = Plot[f[z] - spline /. sol[[1]], {z, a, b}];

El caso cuadrático (clase (C^1))

En esta sección se muestra cómo obtener un spline cuadrático de interpolación una vez especificados los nodos de interpolación, (Subscript[x, 1]) < (Subscript[x, 2]) < ... < (Subscript[x, n]), y los valores que deben ser interpolados,  (Subscript[y, 1]),  (Subscript[y, 2]),  ...,  (Subscript[y, n]).

Si se usan potencias truncadas cuadráticas, el spline buscado se representa como sigue:

spline = 
      (∑_(i = 0)^2) c[i] (z^i) + 
        (∑_(i = 3)^n) c[i] potenctruncada[z - x[[i - 1]], 2];

Imponemos las condiciones de interpolación:

sistema = 
      Table[(spline /. z -> x[[i]]) == y[[i]], {i, n}];

Sabemos que se necesita un dato adicional, usalmente se considera la derivada en uno de los nodos, por ejemplo, el primero. Imponemos el cumplimiento de tal condición y la añadimos al sistema obtenido imponiendo la interpolación en los nodos.

?D

AppendTo[sistema, 
      (D[spline, z] /. z -> x[[1]]) == f'[x[[1]]]]

Definimos las incógnitas del problema, y resolvemos el sistema obtenido:

\(incognitas = Table[c[i], {i, 0, n}]; 


    sol = Solve[sistema, incognitas]\)

Podemos representar gráficamente el spline resultante

grafspl2 = Plot[spline /. sol[[1]], {z, a, b}];

y mostrar conjuntamente su gráfica y la de la función interpolada.

Show[graf, grafspl2];

Asimismo, representamos el error de interpolación.

graferror = Plot[f[z] - spline /. sol[[1]], {z, a, b}];

El caso cúbico (clase (C^2))

Mantenemos el ejemplo de las secciones anteriores y calcularemos diferentes splines cúbicos de clase dos utilizando potencias truncadas cúbicas.

El spline se escibe como combinación lineal de la base elegida del espacio de funciones spline cúbicas de clase dos.

spline = 
      (∑_(i = 0)^3) c[i] (z^i) + 
        (∑_(i = 4)^(n + 1))c[i] 
            potenctruncada[z - x[[i - 2]], 3];

Y se determinan las ecuaciones lineales que garantizan la interpolación de los datos.

sistema = 
      Table[(spline /. z -> x[[i]]) == y[[i]], {i, n}];

A este sistema hay que añadirle dos condiciones adicionales, que dependerán del tipo de spline deseado.

Caso cúbico natural

Añadimos a las ecuaciones anteriores las que se obtiene tras imponer que la derivada segunda en los extremos sea nula.

sistemanatural = 
      Append[sistema, {(D[spline, {z, 2}] /. z -> x[[1]]) == 0, 
            (D[spline, {z, 2}] /. z -> x[[n]]) == 0}] // Flatten

Definimos la lista de incógnitas y resolvemos el sistema.

\(incognitas = Table[c[i], {i, 0, n + 1}]; 

    sol = Solve[sistemanatural, incognitas]\)

Representamos el spline natural calculado y, posteriormente, mostramos su gráfica y la de la función interpolada, y el error de interpolación.

grafsplnatural = Plot[spline /. sol[[1]], {z, a, b}];

Show[graf, grafsplnatural];

graferror = Plot[f[z] - spline /. sol[[1]], {z, a, b}];

Caso cúbico sujeto

El caso cúbico sujeto se trata igualmente. Damos el código necesario.

\(d1 = \(f'\)[x[\([1]\)]]; \ dn = \(f'\)[x[\([n]\)]]; 


    sistemasujeto = 
      Append[sistema, {\((D[spline, z] /. z -> x[\([1]\)])\) == d1, 
            \((D[spline, z] /. z -> x[\([n]\)])\) == dn}] // Flatten
        
\)\(incognitas = Table[c[i], {i, 0, n + 1}]; 


    sol = Solve[sistemasujeto, incognitas]\)

Representamos el spline, la función interpolada y el error de interpolación.

grafspl3sujeto = Plot[spline /. sol[[1]], {z, a, b}];

Show[graf, grafspl3sujeto];

graferror = Plot[f[z] - spline /. sol[[1]], {z, a, b}];

Caso cúbico periódico

Finalizamos considerando el caso periódico.

\(sistemaperiodico = 
      Append[sistema, {
            \((D[spline, z] /. z -> x[\([1]\)])\) == 
              \((D[spline, z] /. z -> x[\([n]\)])\), 
            \((D[spline, {z, 2}] /. z -> x[\([1]\)])\) == 
              \((D[spline, {z, 2}] /. z -> x[\([n]\)])\)}] // Flatten; 


    incognitas = Table[c[i], {i, 0, n + 1}]; 


    sol = Solve[sistemaperiodico, incognitas]\)

grafsplperiodico = Plot[spline /. sol[[1]], {z, a, b}];

Show[graf, grafsplperiodico];

graferror = Plot[f[z] - spline /. sol[[1]], {z, a, b}];

Interpolación spline cúbica de Hermite de orden uno (y por tanto clase (C^1))

Una situación diferente se tiene si hay que construir una función spline cúbica conocidos sus valores y los de su derivada primera en los nodos de interpolación, pues aparece un problema de interpolación de Hermite (de primer orden). La construcción del spline es, en este caso, muy simple, pues se puede determinar intervalo a intervalo resolviendo un problema de interpolación de Hermite polinómico. El interpolante así definido es sólo de clase (C^1), no (C^2).

Desarrollamos un ejemplo, manteniendo los nodos de interpolación y la función que ya hemos empleado anteriormente.

Sólo temenos que formar una lista con los valores de la derivada de la función en los nodos.

d = f'[x]

Representamos los datos que se interpolan asociando a cada punto un segmento de recta que representa a la tangente en el mismo.

\(puntos = 
      ListPlot[Table[{x[[k]], 
            y[[k]]}, {k, n}], 
        PlotStyle -> {PointSize[0.02]}, DisplayFunction -> Identity]; 

    longmin = Min[h]; 

    tangentes = 
      Table[Plot[
          y[[k]] + 
            d[[k]] 
              \((t - x[[k]])\), {t, 
            x[[k]] - \(1\/4\) longmin, 
            x[[k]] + \(1\/4\) longmin}, 
          \ PlotStyle -> {RGBColor[0, 0, 1]}, DisplayFunction -> Identity], {
          k, n}]; 


    Show[Flatten[{puntos, tangentes}], DisplayFunction -> $DisplayFunction]; 
    \)

Se ha usado la orden Flatten[{puntos,tangentes}], que tiene como efecto reducir a una lista unidimensional lo que es una compuesta por una lista de n elementos -tangentes- y un elemento adicional que no forma parte de ella -puntos-. (Nótese que el resultado es una lista de gráficas).

Una vez representada la situación, procedemos a calcular la función spline.

Hallamos las expresiones de los segmentos del spline, haciendo uso de la orden InterpolatingPolynomial.

splH = 
      Table[Expand[
          InterpolatingPolynomial[{{
                x[[k]], {
                  y[[k]], 
                  d[[k]]}}, {
                x[[k + 1]], {
                  y[[k + 1]], 
                  d[[k + 1]]}}}, t]], {
          k, n - 1}]

Generamos sus gráficas individualmente, sin mostrarlas, y con ellas formamos una lista, grafsplH.

\(grafsplH = {}; 


    For[k = 1, k ≤ n - 1, \(k++\), 
	\ 
      AppendTo[grafsplH, 
		\ \ \ \ \ \ \ \ \ \ \ \ \ 
        Plot[splH[[k]], {t, 
            x[[k]], 
            x[[k + 1]]}, 
          PlotStyle -> {Dashing[{0.01, 0.01}]}, 
          DisplayFunction -> Identity]]
\ \ \ \ \ \ \ ]\)

Mostramos la gráfica del interpolante global.

Show[grafsplH, DisplayFunction -> $DisplayFunction];

La mostramos conjuntamente con los datos, las tangentes y la gráfica de la función interpolada.

Show[{grafsplH, puntos, tangentes, graf}, 
      DisplayFunction -> $DisplayFunction];

Actuamos como en la interpolación spline de Lagrange para apreciar el error cometido.

\(errorH = 
      Table[f[t] - splH[[k]], {k, 
          n - 1}]; 

graferrorH = {}; 


    For[k = 1, k ≤ n - 1, \(k++\), 
	\ 
      AppendTo[graferrorH, 
		\ \ \ \ \ \ \ \ \ \ \ \ \ 
        Plot[errorH[[k]], {t, 
            x[[k]], 
            x[[k + 1]]}, 
          PlotStyle -> {Dashing[{0.01, 0.01}]}, 
          DisplayFunction -> Identity]]
\ \ \ \ \ \ \ ]
\)Show[graferrorH, DisplayFunction -> $DisplayFunction];

Ejercicios

1.- Halle el spline cuadrático lagrangiano que interpola la función f(x) = (e^x) sen x en los puntos (Subscript[x, i]) = (π/6) i, i = 0, ..., 6.

2.- Halle el spline cúbico natural que interpola la función f(t) = (((6/5))/((t^2) + 1 + (3/5))) + (1/5)  en los nodos resultantes al dividir el intervalo [-(39/8), (39/8)] en siete partes iguales.

3.- En las mismas condiciones del ejercicio anterior, halle el spline cúbico de interpolación de Hermite de primer orden.

4.- Represente los errores de interpolación correspondientes a los tres problemas anteriores.