8: Interpolación

Considera el siguiente problema: Dados los valores de una función conocida\(y=f(x)\) en una secuencia de puntos ordenados\(x_{0}, x_{1}, \ldots, x_{n}\), encuentra\(f(x)\) para arbitrario\(x .\) Cuando\(x_{0} \leq x \leq x_{n}\), el problema se llama interpolación. Cuando\(x<x_{0}\) o\(x>x_{n}\), el problema se llama extrapolación.

Con\(y_{i}=f\left(x_{i}\right)\), el problema de la interpolación es básicamente uno de dibujar una curva suave a través de los puntos conocidos\(\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots,\left(x_{n}, y_{n}\right) .\) Este no es el mismo problema que dibujar una curva suave que se aproxime a un conjunto de puntos de datos que tienen error experimental. Este último problema se denomina aproximación de mínimos cuadrados, que se considera en el siguiente capítulo.

Es posible interpolar los puntos\(n+1\) conocidos por un polinomio único de grado\(n\). Cuando\(n=1\), el polinomio es una función lineal; cuando\(n=2\), el polinomio es una función cuadrática. Aunque a veces se utilizan polinomios de orden bajo cuando el número de puntos es escaso, los interpolados polinomios de orden superior tienden a ser inestables y no son de mucha utilidad práctica.

Aquí, consideraremos la interpolación polinómica por piezas más útil. Los dos más utilizados son la interpolación lineal por tramos y la interpolación spline cúbica. El primero hace uso de polinomios lineales, y los segundos polinomios cúbicos.

Interpolación lineal por tramos

Aquí, utilizamos polinomios lineales. Esta es la interpolación predeterminada que normalmente se usa al trazar datos.

Supongamos que la función de interpolación es\(y=g(x)\), y como antes, hay\(n+1\) puntos para interpolar. Construimos la función\(g(x)\) a partir de polinomios lineales\(n\) locales. Escribimos

\[g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1}, \nonumber \]

donde

\[g_{i}(x)=a_{i}\left(x-x_{i}\right)+b_{i} \nonumber \]

y\(i=0,1, \ldots, n-1\)

Ahora requerimos\(y=g_{i}(x)\) pasar por los puntos finales\(\left(x_{i}, y_{i}\right)\) y\(\left(x_{i+1}, y_{i+1}\right) .\) tenemos

\[\begin{aligned} y_{i} &=b_{i} \\ y_{i+1} &=a_{i}\left(x_{i+1}-x_{i}\right)+b_{i} \end{aligned} \nonumber \]

Por lo tanto,\(g_{i}(x)\) se determina que los coeficientes de

\[a_{i}=\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}, \quad b_{i}=y_{i} \nonumber \]

Aunque la interpolación lineal por tramos es ampliamente utilizada, particularmente en rutinas de trazado, sufre de una discontinuidad en la derivada en cada punto. Esto da como resultado una función que puede no parecer suave si los puntos están demasiado espaciados. A continuación consideramos un algoritmo más desafiante que utiliza polinomios cúbicos.

Interpolación spline cúbica

Los\(n+1\) puntos a ser interpolados son de nuevo

\[\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots\left(x_{n}, y_{n}\right) . \nonumber \]

Aquí, utilizamos\(n\) polinomios cúbicos por partes para la interpolación,

\[g_{i}(x)=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i}, \quad i=0,1, \ldots, n-1, \nonumber \]

con la función de interpolación global escrita como

\[g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1} . \nonumber \]

Para lograr una interpolación suave lo imponemos\(g(x)\) y sus derivadas primera y segunda son continuas. El requisito que\(g(x)\) es continuo (y pasa por todos los\(n+1\) puntos) da como resultado las dos restricciones

\[\begin{align} g_{i}\left(x_{i}\right) &=y_{i}, \quad i=0 \text { to } n-1 \\ g_{i}\left(x_{i+1}\right) &=y_{i+1}, \quad i=0 \text { to } n-1 \end{align} \nonumber \]

El requisito que\(g^{\prime}(x)\) es continuo da como resultado

\[g_{i}^{\prime}\left(x_{i+1}\right)=g_{i+1}^{\prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 \nonumber \]

Y el requerimiento que\(g^{\prime \prime}(x)\) es continuo da como resultado

\[g_{i}^{\prime \prime}\left(x_{i+1}\right)=g_{i+1}^{\prime \prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 . \nonumber \]

Hay polinomios\(n\) cúbicos\(g_{i}(x)\) y cada polinomio cúbico tiene cuatro coeficientes libres; por lo tanto, hay un total de coeficientes\(4 n\) desconocidos. El número de ecuaciones de restricción de (8.7) - (8.10) es\(2 n+2(n-1)=4 n-2\). Con\(4 n-2\) restricciones e\(4 n\) incógnitas, se requieren dos condiciones más para una solución única. Estas suelen elegirse para ser condiciones extra en el primer\(g_{0}(x)\) y último\(g_{n-1}(x)\) polinomios. Discutiremos estas condiciones adicionales más adelante.

Ahora se procede a determinar ecuaciones para los coeficientes desconocidos de los polinomios cúbicos. Los polinomios y sus dos primeras derivadas están dados por

\[\begin{align} g_{i}(x) &=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i} \\ g_{i}^{\prime}(x) &=3 a_{i}\left(x-x_{i}\right)^{2}+2 b_{i}\left(x-x_{i}\right)+c_{i} \\ g_{i}^{\prime \prime}(x) &=6 a_{i}\left(x-x_{i}\right)+2 b_{i} \end{align} \nonumber \]

Consideraremos a su vez las cuatro condiciones (8.7) - (8.10). Desde (8.7) y (8.11), tenemos

\[d_{i}=y_{i}, \quad i=0 \text { to } n-1, \nonumber \]

que resuelve directamente para todos los\(d\) coeficientes.

Para satisfacer (8.8), primero definimos

\[h_{i}=x_{i+1}-x_{i}, \nonumber \]

y

\[f_{i}=y_{i+1}-y_{i} . \nonumber \]

Ahora, a partir de (8.8) y (8.11), usando (8.14), obtenemos las\(n\) ecuaciones

\[a_{i} h_{i}^{3}+b_{i} h_{i}^{2}+c_{i} h_{i}=f_{i}, \quad i=0 \text { to } n-1 . \nonumber \]

De (8.9) y (8.12) obtenemos las\(n-1\) ecuaciones

\[3 a_{i} h_{i}^{2}+2 b_{i} h_{i}+c_{i}=c_{i+1}, \quad i=0 \text { to } n-2 \nonumber \]

De (8.10) y (8.13) obtenemos las\(n-1\) ecuaciones

\[3 a_{i} h_{i}+b_{i}=b_{i+1} \quad i=0 \text { to } n-2 \text {. } \nonumber \]

Será útil incluir una definición del coeficiente\(b_{n}\), que ahora falta. (El índice de los coeficientes polinomiales cúbicos sólo va hasta\(n-1\).) Simplemente extendemos (8.19) hasta\(i=n-1\) y así escribimos

\[3 a_{n-1} h_{n-1}+b_{n-1}=b_{n}, \nonumber \]

que puede ser visto como una definición de\(b_{n}\).

Ahora se procede a eliminar los conjuntos de coeficientes\(a\) - y c (con los\(d\) -coeficientes ya eliminados en (8.14)) para encontrar un sistema de ecuaciones lineales para los\(b\) -coeficientes. A partir de (8.19) y\((8.20)\), podemos resolver para los\(n a\) -coeficientes para encontrar

\[a_{i}=\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right), \quad i=0 \text { to } n-1 . \nonumber \]

A partir de (8.17), podemos resolver para los\(n\) coeficientes c de la siguiente manera:

\[\begin{align} \nonumber c_{i} &=\frac{1}{h_{i}}\left(f_{i}-a_{i} h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ \nonumber &=\frac{1}{h_{i}}\left(f_{i}-\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right) h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ &=\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right), \quad i=0 \text { to } n-1 \end{align} \nonumber \]

Ahora podemos encontrar una ecuación para los\(b\) -coeficientes sustituyendo\((8.21)\) y\((8.22)\) en\((8.18)\):

\[\begin{array}{r} 3\left(\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right)\right) h_{i}^{2}+2 b_{i} h_{i}+\left(\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right)\right) \\ =\left(\frac{f_{i+1}}{h_{i+1}}-\frac{1}{3} h_{i+1}\left(b_{i+2}+2 b_{i+1}\right)\right), \nonumber \end{array} \nonumber \]

lo que simplifica a

\[\frac{1}{3} h_{i} b_{i}+\frac{2}{3}\left(h_{i}+h_{i+1}\right) b_{i+1}+\frac{1}{3} h_{i+1} b_{i+2}=\frac{f_{i+1}}{h_{i+1}}-\frac{f_{i}}{h_{i}} \nonumber \]

una ecuación que es válida para\(i=0\) a\(n-2\). Por lo tanto, (8.23) representan\(n-1\) ecuaciones para los\(b\) coeficientes\(n+1\) desconocidos. En consecuencia, escribimos la ecuación matricial para los\(b\) -coeficientes, dejando ausentes la primera y última fila, como

\[\left(\begin{array}{cccccccc} \ldots & \ldots & \ldots & \ldots & \text { missing } & \ldots & \ldots \\ \frac{1}{3} h_{0} & \frac{2}{3}\left(h_{0}+h_{1}\right) & \frac{1}{3} h_{1} & \ldots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & \frac{1}{3} h_{n-2} & \frac{2}{3}\left(h_{n-2}+h_{n-1}\right) & \frac{1}{3} h_{n-1} \\ \ldots & \ldots & \ldots & \ldots & \text { missing } & \ldots & \ldots \end{array}\right)\left(\begin{array}{c} b_{0} \\ b_{1} \\ \vdots \\ b_{n-1} \\ b_{n} \end{array}\right)\left(\begin{array}{c} \text { missing } \\ \frac{f_{1}}{h_{1}}-\frac{f_{0}}{h_{0}} \\ \vdots \\ \frac{f_{n-1}}{h_{n-1}}-\frac{f_{n-2}}{h_{n-2}} \\ \text { missing } \end{array}\right) . \nonumber \]

Una vez que se especifican las ecuaciones primera y última faltantes, la ecuación matricial para los\(b\) coeficientes se puede resolver mediante eliminación gaussiana. Y una vez determinados los\(b\) -coeficientes, los\(c\) coeficientes\(a\) - y -también se pueden determinar a partir de (8.21) y (8.22), con los\(d\) -coeficientes ya conocidos a partir de (8.14). Los polinomios cúbicos por tramos, entonces, son conocidos y\(g(x)\) pueden ser utilizados para la interpolación a cualquier valor que\(x\) satisfaga\(x_{0} \leq x \leq x_{n}\).

La primera y última ecuaciones faltantes se pueden especificar de varias maneras, y aquí mostramos las dos formas que están permitidas por la función de MATLAB spline.m La primera forma se debe usar cuando la derivada\(g^{\prime}(x)\) es conocida en los puntos finales\(x_{0}\) y es\(x_{n} ;\) decir, supongamos que conocemos los valores de\(\alpha\) y\(\beta\) tal

\[g_{0}^{\prime}\left(x_{0}\right)=\alpha, \quad g_{n-1}^{\prime}\left(x_{n}\right)=\beta . \nonumber \]

A partir del valor conocido de\(\alpha\), y usando (8.12) y (8.22), tenemos

\[\begin{aligned} \alpha &=c_{0} \\ &=\frac{f_{0}}{h_{0}}-\frac{1}{3} h_{0}\left(b_{1}+2 b_{0}\right) \end{aligned} \nonumber \]

Por lo tanto, se determina que la primera ecuación que falta es

\[\frac{2}{3} h_{0} b_{0}+\frac{1}{3} h_{0} b_{1}=\frac{f_{0}}{h_{0}}-\alpha \nonumber \]

Desde el valor conocido de\(\beta\), y usando\((8.12),(8.21)\), y\((8.22)\), tenemos

\[\begin{aligned} \beta &=3 a_{n-1} h_{n-1}^{2}+2 b_{n-1} h_{n-1}+c_{n-1} \\ &=3\left(\frac{1}{3 h_{n-1}}\left(b_{n}-b_{n-1}\right)\right) h_{n-1}^{2}+2 b_{n-1} h_{n-1}+\left(\frac{f_{n-1}}{h_{n-1}}-\frac{1}{3} h_{n-1}\left(b_{n}+2 b_{n-1}\right)\right), \end{aligned} \nonumber \]

lo que simplifica a

\[\frac{1}{3} h_{n-1} b_{n-1}+\frac{2}{3} h_{n-1} b_{n}=\beta-\frac{f_{n-1}}{h_{n-1}} \nonumber \]

para ser utilizada como la última ecuación faltante.

La segunda forma de especificar la primera y última ecuaciones faltantes se llama la condición no-un nudo, que asume que

\[g_{0}(x)=g_{1}(x), \quad g_{n-2}(x)=g_{n-1}(x) \nonumber \]

Considerando la primera de estas ecuaciones, desde (8.11) tenemos

\[a_{0}\left(x-x_{0}\right)^{3}+b_{0}\left(x-x_{0}\right)^{2}+c_{0}\left(x-x_{0}\right)+d_{0} \nonumber \]

\[=a_{1}\left(x-x_{1}\right)^{3}+b_{1}\left(x-x_{1}\right)^{2}+c_{1}\left(x-x_{1}\right)+d_{1} . \nonumber \]

Ahora se puede probar que dos polinomios cúbicos son idénticos si a algún valor de\(x\), los polinomios y sus tres primeras derivadas son idénticos. Nuestras condiciones de continuidad en\(x=x_{1}\) ya requieren que a este valor de\(x\) estos dos polinomios y sus dos primeras derivadas sean idénticos. Los polinomios en sí serán idénticos, entonces, si sus terceras derivadas también son idénticas en\(x=x_{1}\), o si

\[a_{0}=a_{1} . \nonumber \]

Desde (8.21), tenemos

\[\frac{1}{3 h_{0}}\left(b_{1}-b_{0}\right)=\frac{1}{3 h_{1}}\left(b_{2}-b_{1}\right), \nonumber \]

o después de la simplificación

\[h_{1} b_{0}-\left(h_{0}+h_{1}\right) b_{1}+h_{0} b_{2}=0, \nonumber \]

lo que nos proporciona nuestra primera ecuación faltante. Un argumento similar en\(x=x_{n-1}\) también nos proporciona nuestra última ecuación,

\[h_{n-1} b_{n-2}-\left(h_{n-2}+h_{n-1}\right) b_{n-1}+h_{n-2} b_{n}=0 . \nonumber \]

Las subrutinas de MATLAB spline.m y ppval.m se pueden utilizar para interpolación spline cúbica (ver también interp1.m). Ilustraré estas rutinas en clase y publicaré código de muestra en el sitio web del curso.

Interpolación multidimensional

Supongamos que estamos interpolando el valor de una función de dos variables,

\[z=f(x, y) . \nonumber \]

Los valores conocidos vienen dados por

\[z_{i j}=f\left(x_{i}, y_{j}\right), \nonumber \]

con\(i=0,1, \ldots, n\) y\(j=0,1, \ldots, n .\) Tenga en cuenta que los\((x, y)\) puntos en\(f(x, y)\) los que se conocen se encuentran en una cuadrícula en el\(x-y\) plano.

\(z=g(x, y)\)Sea la función de interpolación, satisfaciendo\(z_{i j}=g\left(x_{i}, y_{j}\right) .\) Una interpolación bidimensional para encontrar el valor de\(g\) en el punto se\((x, y)\) puede hacer primero realizando interpolaciones\(n+1\) unidimensionales en \(y\)para encontrar el valor de\(g\) en los\(n+1\) puntos\(\left(x_{0}, y\right),\left(x_{1}, y\right), \ldots\),\(\left(x_{n}, y\right)\), seguido de una sola interpolación unidimensional en\(x\) para encontrar el valor de\(g\) at\((x, y)\).

En otras palabras, la interpolación bidimensional en una cuadrícula de dimensión\((n+1) \times(n+1)\) se realiza primero realizando interpolaciones\(n+1\) unidimensionales al valor\(y\) seguido de una sola interpolación unidimensional al valor\(x\). La interpolación bidimensional puede generalizarse a dimensiones superiores. Las funciones de MATLAB para realizar interpolación bidimensional y tridimensional son interp2.m e interp3.m.