Interpolación
El problema de interpolación, el condicionamiento de la prescripción que produce el interpolante y la interpolación lineal a trozos.
Pasamos ahora al capítulo 5 del libro de Driscoll y Braun, dedicado a la interpolación. Seguimos las secciones 5.1 y 5.2.
1 El problema de interpolación
Dados \(n+1\) puntos \((t_0,y_0),\dots,(t_n,y_n)\) con los \(t_k\) distintos, llamados nodos, el problema de interpolación consiste en encontrar una función \(p\), llamada interpolante, tal que
\[p(t_k)=y_k,\qquad k=0,\dots,n.\]
Tomemos los tres puntos \((0,1)\), \((1,3)\) y \((2,2)\). Hay muchas funciones que pasan por ellos, y vemos tres.
Lineal a trozos. En cada subintervalo usamos una recta que une sus dos extremos:
\[p_1(x)=1+2x\ \text{ en }[0,1],\qquad p_2(x)=4-x\ \text{ en }[1,2].\]
Cuadrática a trozos. Usamos ahora un polinomio de grado a lo más \(2\) en cada subintervalo,
\[p_1(x)=a_0+a_1x+a_2x^2\ \text{ en }[0,1],\qquad p_2(x)=b_0+b_1x+b_2x^2\ \text{ en }[1,2].\]
Pasar por los cuatro extremos impone cuatro condiciones sobre seis coeficientes, así que quedan dos grados de libertad. Los podemos usar para que el empalme en \(x=1\) sea suave, \(p_1'(1)=p_2'(1)\), y para fijar, por ejemplo, \(p_1'(0)=0\). De \(p_1(0)=1\), \(p_1'(0)=0\) y \(p_1(1)=3\) sale \(p_1(x)=1+2x^2\), con \(p_1'(1)=4\). Escribiendo \(p_2\) alrededor de \(x=1\) con esa misma pendiente, \(p_2(x)=3+4(x-1)+b_2(x-1)^2\), y la condición \(p_2(2)=2\) da \(b_2=-5\):
\[p_1(x)=1+2x^2,\qquad p_2(x)=3+4(x-1)-5(x-1)^2.\]
La condición \(p_1'(0)=0\) es solo una elección posible, otra elección da otra cuadrática a trozos.
Una sola cuadrática. También podemos usar un único polinomio de grado \(2\) para todo el intervalo, el único que pasa por los tres puntos:
\[p(x)=1+\frac72\,x-\frac32\,x^2.\]
En general, \(n+1\) puntos se pueden interpolar con un único polinomio de grado \(n\). Sin embargo, esto puede generar oscilaciones muy grandes en los extremos del intervalo, con errores y sensibilidad altos a los datos (fenómeno de Runge). Por eso es preferible interpolar con polinomios de grado pequeño en cada subintervalo. Los ejemplos de la sección 5.1 del libro lo muestran.
2 Condicionamiento de la interpolación
Fijados los nodos, una prescripción de interpolación es una función \(I\) que a cada vector de datos \(y=(y_0,\dots,y_n)\) le asigna su interpolante \(p=I(y)\), es decir, \(p(t_k)=y_k\). Solo consideraremos prescripciones lineales:
\[I(\alpha y+\beta z)=\alpha\,I(y)+\beta\,I(z)\]
para todo par de vectores \(y,z\) y todo par de escalares \(\alpha,\beta\).
Con un ejemplo se ve más claro. Sean \(t_0=0\), \(t_1=1\), \(t_2=2\), e \(I\) la prescripción que a los datos \((0,y_0)\), \((1,y_1)\), \((2,y_2)\) les asigna la función que une \((0,y_0)\) con \((1,y_1)\) mediante una recta y \((1,y_1)\) con \((2,y_2)\) mediante una cuadrática, de modo que en total \(p\) sea de clase \(C^1\). La recta tiene pendiente \(y_1-y_0\), así que la cuadrática debe salir de \((1,y_1)\) con esa pendiente y pasar por \((2,y_2)\). El resultado es
\[p(x)=I(y)(x)=\begin{cases} y_0+(y_1-y_0)\,x, & 0\le x\le 1,\\[2mm] y_1+(y_1-y_0)(x-1)+(y_0-2y_1+y_2)(x-1)^2, & 1\le x\le 2. \end{cases}\]
¿Qué pasa si calculamos \(I(\alpha y+\beta z)(x)\)? Para \(x\in[0,1]\),
\[\begin{aligned} I(\alpha y+\beta z)(x)&=(\alpha y_0+\beta z_0)+\big[(\alpha y_1+\beta z_1)-(\alpha y_0+\beta z_0)\big]\,x\\ &=\alpha\big[y_0+(y_1-y_0)x\big]+\beta\big[z_0+(z_1-z_0)x\big]\\ &=\alpha\,I(y)(x)+\beta\,I(z)(x), \end{aligned}\]
y de manera similar si \(x\in[1,2]\).
Como \(I\) es lineal en \(y\), podemos escribir el vector de datos como combinación lineal de los vectores canónicos, \(y=\sum_{k=0}^n y_k\,e_k\), donde \(e_k\) es la columna \(k\) de la matriz identidad (numerada desde \(0\)). Por linealidad,
\[I(y)=\sum_{k=0}^{n}y_k\,I(e_k).\]
Las funciones \(\varphi_k=I(e_k)\) son las que se llaman funciones cardinales: \(\varphi_k\) es el interpolante de los datos que valen \(1\) en el nodo \(t_k\) y \(0\) en los demás, es decir, \(\varphi_k(t_j)=1\) si \(j=k\) y \(0\) si \(j\neq k\). En el ejemplo, evaluando la fórmula de \(I\) en \(e_0=(1,0,0)\), \(e_1=(0,1,0)\) y \(e_2=(0,0,1)\) se obtiene
\[I(e_0)(x)=\begin{cases}1-x, & x\in[0,1],\\ -(x-1)+(x-1)^2, & x\in[1,2],\end{cases}\qquad I(e_1)(x)=\begin{cases}x, & x\in[0,1],\\ 1+(x-1)-2(x-1)^2, & x\in[1,2],\end{cases}\]
\[I(e_2)(x)=\begin{cases}0, & x\in[0,1],\\ (x-1)^2, & x\in[1,2],\end{cases}\]
y por lo tanto \(I(y)=y_0\,I(e_0)+y_1\,I(e_1)+y_2\,I(e_2)\).
Medimos los datos con la norma \(\|d\|_\infty=\max_k|d_k|\) y las funciones con \(\|p\|_\infty=\max_{t_0\le x\le t_n}|p(x)|\).
El número de condicionamiento absoluto de \(I\) es
\[\kappa(I)=\max_{d\neq0}\frac{\|I(y+d)-I(y)\|_\infty}{\|d\|_\infty},\]
es decir, el máximo del cociente entre el error (absoluto) en el resultado y el error (absoluto) en el dato.
El número de condicionamiento absoluto de \(I\) satisface
\[\max_{k}\|\varphi_k\|_\infty\;\le\;\kappa(I)\;\le\;\sum_{k=0}^{n}\|\varphi_k\|_\infty.\tag{1}\]
Demostración. Por linealidad, \(I(y+d)-I(y)=I(d)\), y entonces
\[\frac{\|I(y+d)-I(y)\|_\infty}{\|d\|_\infty}=\frac{\|I(d)\|_\infty}{\|d\|_\infty}=\left\|I\!\left(\frac{d}{\|d\|_\infty}\right)\right\|_\infty .\]
El número de condicionamiento es el máximo de esta cantidad sobre todos los \(d\neq0\) posibles, y no depende de \(y\). Sea \(k^*\) un índice en el que \(\max_k\|\varphi_k\|_\infty\) se alcanza, es decir, \(\|\varphi_{k^*}\|_\infty=\max_k\|\varphi_k\|_\infty\). Tomando \(d=e_{k^*}\) se tiene \(\|d\|_\infty=1\) e \(I(d)=\varphi_{k^*}\), de modo que el cociente vale \(\|\varphi_{k^*}\|_\infty\). Como \(\kappa(I)\) es el máximo sobre todos los \(d\), resulta \(\kappa(I)\ge\|\varphi_{k^*}\|_\infty\), que es la desigualdad de la izquierda en (1). La otra desigualdad es consecuencia directa de la desigualdad triangular y se deja como ejercicio. \(\blacksquare\)
Demuestre la desigualdad de la derecha en (1) y resuelva los ejercicios 5.1.4 y 5.1.5 de Fundamentals of Numerical Computation, disponibles en https://fncbook.com/julia/interpolation/#exercises.
3 Interpolación lineal a trozos
Consideremos ahora los nodos ordenados, \(t_0<t_1<\dots<t_n\), y la interpolación lineal a trozos (PL, por piecewise linear): en cada subintervalo \([t_k,t_{k+1}]\) el interpolante es el segmento que une \((t_k,y_k)\) con \((t_{k+1},y_{k+1})\),
\[p(x)=y_k+\frac{y_{k+1}-y_k}{t_{k+1}-t_k}\,(x-t_k),\qquad x\in[t_k,t_{k+1}],\quad k=0,\dots,n-1.\]
Las funciones cardinales de PL se llaman funciones sombrero, por su forma:
\[H_k(x)=\begin{cases} \dfrac{x-t_{k-1}}{t_k-t_{k-1}}, & x\in[t_{k-1},t_k],\\[3mm] \dfrac{t_{k+1}-x}{t_{k+1}-t_k}, & x\in[t_k,t_{k+1}],\\[3mm] 0, & \text{en otro caso.} \end{cases}\]
Para \(H_0\) no hay rama izquierda y para \(H_n\) no hay rama derecha, porque no existen \(t_{-1}\) ni \(t_{n+1}\).
Estas funciones son continuas en \(\mathbb{R}\) y lineales a trozos, y son cardinales para PL: \(H_k(t_j)=1\) si \(j=k\) y \(0\) si \(j\neq k\). Entonces
\[p(x)=\sum_{k=0}^{n}y_k\,H_k(x).\]
El siguiente tablero dibuja los sombreros, los sombreros multiplicados por los datos y el interpolante que resulta de sumarlos. Se pueden arrastrar los puntos para cambiar los datos, cambiar el número de subintervalos y perturbar los nodos para que no sean equiespaciados, y leer en una tabla el valor de cada sombrero en un punto.
Demuestre que el número de condicionamiento absoluto de la interpolación lineal a trozos es \(\kappa(\mathrm{PL})=1\).
3.1 Por qué condicionamiento absoluto
El adjetivo absoluto significa que el error en el dato y el error en el resultado se miden tal cual, sin dividir por el tamaño de \(y\) ni por el de \(I(y)\). Entonces \(\kappa(I)\) dice cuánto se amplifica, a lo sumo, un error absoluto: si cada \(y_k\) cambia en no más de \(\delta\), el interpolante cambia, en la norma del máximo, en no más de \(\kappa(I)\,\delta\). Lo usamos por dos razones.
La primera es que, por ser \(I\) lineal, la diferencia \(I(y+d)-I(y)=I(d)\) no depende de \(y\), así que \(\kappa(I)\) es un único número que depende solo de los nodos y del método de interpolación, y vale como garantía para todos los datos a la vez. La segunda es que el condicionamiento relativo en un dato \(y\neq0\) sí depende de \(y\),
\[\kappa_{\text{rel}}(y)=\kappa(I)\,\frac{\|y\|_\infty}{\|I(y)\|_\infty},\]
pero como el interpolante pasa por los nodos se tiene \(\|I(y)\|_\infty\ge\|y\|_\infty\), y por lo tanto \(\kappa_{\text{rel}}(y)\le\kappa(I)\) para todo \(y\). Acotar el condicionamiento absoluto controla también el relativo.
4 Convergencia de la interpolación lineal a trozos
Sea \(f\) con segunda derivada continua en \([a,b]\), sea \(h=\dfrac{b-a}{n}\) y sea \(p_n\) el interpolante PL de los datos \((t_i,f(t_i))\), con nodos equiespaciados \(t_i=a+ih\). Entonces
\[\|f-p_n\|_\infty\le\|f''\|_\infty\,h^2.\]
Como \(h=(b-a)/n\), esto dice que la convergencia es algebraica, \(\|f-p_n\|_\infty=O(h^2)=O(n^{-2})\), y si el error no es \(O(h^3)\), entonces \(2\) es el orden de precisión.
4.1 El orden aquí y en los métodos iterativos
La palabra orden aparece en los dos contextos, pero mide cosas distintas. En los métodos iterativos (punto fijo, Newton, secante) hay una sola sucesión de errores \(\varepsilon_k=x_k-x^*\), uno por cada iteración, y el orden \(\alpha\) relaciona cada error con el anterior:
\[\lim_{k\to\infty}\frac{|\varepsilon_{k+1}|}{|\varepsilon_k|^\alpha}=L>0.\]
Con \(\alpha=1\) la convergencia es lineal, con \(\alpha=2\) es cuadrática. Lo que va cambiando es el número de pasos \(k\), y el error baja de forma geométrica o más rápido: la convergencia lineal es una recta en una gráfica con el error en escala logarítmica contra \(k\).
En la interpolación no hay iteración ni error anterior. Para cada \(n\) hay un interpolante distinto \(p_n\), y el error \(E(h)=\|f-p_n\|_\infty\) depende del tamaño de la malla. El orden \(q\) es el exponente de \(h\):
\[E(h)\approx C\,h^{q}.\]
Para estimarlo se comparan dos mallas, \(E(h_1)/E(h_2)\approx(h_1/h_2)^q\), y al duplicar \(n\) queda \(q\approx\log_2\big(E(h)/E(h/2)\big)\). Equivale a la pendiente de la recta que se ve al graficar \(E\) contra \(h\) con ambos ejes en escala logarítmica. Para PL, \(q=2\): al dividir \(h\) entre \(2\), el error se divide aproximadamente entre \(4\). Con \(f(x)=\operatorname{sen}x\) en \([0,\pi]\):
| \(n\) | \(h\) | \(\lVert f-p_n\rVert_\infty\) | \(E(2h)/E(h)\) | \(\log_2\) del cociente |
|---|---|---|---|---|
| \(4\) | \(0{,}78540\) | \(7{,}0378\times10^{-2}\) | ||
| \(8\) | \(0{,}39270\) | \(1{,}8846\times10^{-2}\) | \(3{,}734\) | \(1{,}901\) |
| \(16\) | \(0{,}19635\) | \(4{,}7921\times10^{-3}\) | \(3{,}933\) | \(1{,}976\) |
| \(32\) | \(0{,}09817\) | \(1{,}2031\times10^{-3}\) | \(3{,}983\) | \(1{,}994\) |
| \(64\) | \(0{,}04909\) | \(3{,}0109\times10^{-4}\) | \(3{,}996\) | \(1{,}998\) |
| \(128\) | \(0{,}02454\) | \(7{,}5292\times10^{-5}\) | \(3{,}999\) | \(2{,}000\) |
El cociente de errores tiende a \(4\) y el logaritmo en base \(2\) del cociente tiende a \(2\), que es el orden de precisión.
Conviene notar la diferencia de velocidad. Un orden \(2\) de este tipo es mucho más lento que la convergencia cuadrática de Newton, que duplica las cifras correctas en cada paso. Aquí cada vez que duplicamos \(n\), y con ello el costo, el error se divide entre \(4\), es decir, ganamos solo unas \(0{,}6\) cifras decimales. El error decae como una potencia de \(n\) y no como una geométrica en \(k\), y de ahí el nombre de convergencia algebraica.
Resuelva los ejercicios 5.2.2 a 5.2.4 de Fundamentals of Numerical Computation, disponibles en https://fncbook.com/julia/pwlin/#exercises.