Interpolación por el método de Newton

En artículos anteriores he hablado sobre la obtención del polinomio interpolador de grado $n$ ( dado un conjunto de $n+1$ puntos ) mediante el método de determinación de los coeficientes del polinomio pedido, resolviendo el sistema de ecuaciones lineales que aparece al imponer que los datos lo satisfagan, que puedes repasar siguiendo este enlace. También he hablado del método de Lagrange para obtener dicho polinomio, bastante más sofisticado y con claras ventajas de cálculo sobre el primero.

Los métodos de interpolación son importantes por su aplicabilidad numérica a muchos tipos de problemas; por ejemplo a la hora de integrar numéricamente una función cuya primitiva no sea fácil obtener o incluso para los casos en que dicha primitiva no sea una función elemental. Y, a pesar de que la solución obtenida así no será exacta, bien puede controlarse la imprecisión que ello conlleva. La idea clave de la integración numérica consiste en sustituir dicha función por una función aproximada, cuya integración no presente problemas, en un cierto intervalo. Desde luego, estos contenidos se cosideran de ampliación del currículo estándar de bachillerato. Y como ya sabes, suelo incidir en ellos una vez acabado ya el curso. Si tienes intención de cursar un grado en el que las matemáticas tengan un protagonismo relevante, creo sinceramente que estos contenidos de ampliación pueden serte de gran ayuda y de mucha motivación.

En este artículo te voy a hablar de otro método básico de interpolación polinómica: el método de Newton. Te darás cuenta de que el ejemplo que he utilizado en todos los artículos dedicados a explicar la interpolación polinómica ( el m. de determinación de los coeficientes por resolución del sistema de ecuaciones, el de Lagrange, y, ahora, el de Newton ) es el mismo ( con sólo tres puntos como datos, dando lugar a un polinomio de grado 2 ), pues he intentado que se entendieran bien los pasos de cálculo, a partir de un ejemplo sencillo. Sin embargo, te animo a que te plantees otros ejemplos con más de tres puntos ( 5 o 6 puntos para empezar estaría bien ), y que utilices para ello los métodos de Lagrange y de Newton, para evitar ( si utilizas el primer método ) sistemas de ecuaciones lineales de gran dimensión, que, en principio, suponen una gran cantidad de operaciones aritméticas a realizar.

El no tener que recalcularlo todo al añadir más puntos al conjunto de datos es otra razón importante a la hora de recomendarte el método de Newton frente, por ejemplo, al primer método ( el de la resolución del sistema de ecuaciones lineales ), siempre y cuando mantengas la misma distancia entre las abscisas de los nuevos puntos añadidos. Si ésto no fuese posible hacerlo por alguna razón, podrías utilizar entonces el método de Lagrange.

-oOo-
Sin más dilación, entro ya en materia sobre el método de interpolación de Newton. Consideremos el siguiente conjunto de $n+1$ puntos del plano cartesiano:

de modo de las abscisas sean equidistantes, esto es, $x_1-x_0=x_2-x_1=x_3-x_2=\ldots=h$


Llamamos diferencias primeras a las $n$ siguientes:
  $\Delta\,f(x_0)=f(x_1)-f(x_0)$
  $\Delta\,f(x_1)=f(x_2)-f(x_1)$
  $\Delta\,f(x_2)=f(x_3)-f(x_2)$
  $\ldots$
  $\Delta\,f(x_{n-1})=f(x_{n})-f(x_{n-1})$

Denominamos diferencias segundas a las $n-1$ siguientes:
  $\Delta^{2}\,f(x_0)=\Delta\,f(x_1)-\Delta\,f(x_0)$
  $\Delta^{2}\,f(x_1)=\Delta\,f(x_2)-\Delta\,f(x_1)$
  $\Delta^{2}\,f(x_2)=\Delta\,f(x_3)-\Delta\,f(x_2)$
  $\ldots$
  $\Delta^{2}\,f(x_{n-2})=\Delta\,f(x_{n-1})-\Delta\,f(x_{n-2})$

Denominamos diferencias terceras a las $n-2$ siguientes:
  $\Delta^{3}\,f(x_0)=\Delta^{2}\,f(x_1)-\Delta^{2}\,f(x_0)$
  $\Delta^{3}\,f(x_1)=\Delta^{2}\,f(x_2)-\Delta^{2}\,f(x_1)$
  $\Delta^{3}\,f(x_2)=\Delta^{2}\,f(x_3)-\Delta^{2}\,f(x_2)$
  $\ldots$
  $\Delta^{3}\,f(x_{n-3})=\Delta^{2}\,f(x_{n-2})-\Delta^{2}\,f(x_{n-3})$

...

y de manera génerica, diferencias $p$-ésimas ( $1 \le p \prec n$ ) a las $n-p$ restas $$\Delta^{p}\,f(x_{x_i})=\Delta^{p-1}\,f(x_{i+1})-\Delta^{p-1}\,f(x_{i})$$

Vamos a probar ahora que si conocemos $\{f(x_0),\Delta\,f(x_0),\Delta^{2}\,f(x_0),\ldots,\Delta^{k}\,f(x_0)\}$ podemos conocer el valor de $f(x_k)$

En efecto, de las diferencias primeras, podemos escribir $$f(x_1)=f(x_0)+\Delta\,f(x_0)$$ y $$f(x_2)=f(x_1)+\Delta\,f(x_1)=f(x_0)+\Delta\,f(x_0)+\Delta\,f(x_1)$$ y como, por otra parte, de las diferencias segundas, $\Delta\,f(x_1)=\Delta^{2}\,f(x_0)+\Delta\,f(x_0)$, podemos acabar escribiendo que $$f(x_2)=f(x_0)+\Delta^{2}\,f(x_0)+2\,\Delta\,f(x_0)$$

De manera análoga, podremos escribir $f(x_3)$ en términos de $f(x_0)$ y de las diferencias de orden tres dos y uno: $$f(x_3)=f(x_0)+3\,\Delta^{2}\,f(x_0)+3\,\Delta\,f(x_0)+\Delta^{3}\,f(x_0)$$ y, en general ( emerge un patrón claro ): $$\displaystyle f(x_k)=\sum_{i=0}^{k}\,\binom{k}{i}\,\Delta^{i}\,f(x_0)\;\text{para todo}\;k=0,1,2,\ldots,n$$

Vamos a comprobar esta fórmula con el siguiente ejemplo. Sean los siguientes tres puntos del plano

A partir de esta información de partida podemos construir una tabla de diferencias:

que, con los datos del ejemplo, se concreta así:

Como tenemos $3=2+1$ puntos, $n=2$; y $k$ tomará los valores: $0,1$ y $2$. Observemos que reproducimos el interior de la tabla; en efecto:

  $f(x_0)=-1$

  $\displaystyle f(x_1)=\binom{1}{0}\,f(x_0)+\binom{1}{1}\,\Delta\,f(x_0)=1\cdot(-1)+1\cdot (-1)=-2$

  $\displaystyle f(x_2)=\binom{2}{0}\,f(x_0)+\binom{2}{1}\,\Delta\,f(x_0)+ \binom{2}{2}\,\Delta^{2}\,f(x_0)=1\cdot (-1)+2\cdot (-1)+1\cdot 2 = -1$

-oOo-
Pues bien, para $k:=n$ podemos escribir

$$\displaystyle f(x_n)=\sum_{i=0}^{n}\,\binom{n}{i}\,\Delta^{i}\,f(x_0)$$ y desarrollando los coeficientes binomiales podemos expresarlo también así

$$\displaystyle f(x_0)+n\,\Delta\,f(x_0)+\dfrac{n\,(n-1)}{2}\,\Delta^{2}\,f(x_0)+\overset{\underbrace{n+1}}{\ldots}+\dfrac{n\,(n-1)\,(n-2)\,\ldots\,1}{n!}\,\Delta^{n}\,f(x_0)$$

Tengamos en cuenta ahora la equidistancia de las abscisas consecutivas ( recordemos que hemos supuesto que dicha distancia es $h$ ), entonces:

$$x_n=x_0+n\,h \Rightarrow n=\dfrac{x_n-x_0}{h}$$

$$x_n=x_1+(n-1)\,h \Rightarrow n-1=\dfrac{x_n-x_1}{h}$$

$$\ldots$$

$$x_n=x_{n-1}+1\cdot h \Rightarrow 1=\dfrac{x_n-x_{n-1}}{h}$$

así que sustituyendo en la expresión de $f(x_n)$ llegamos a

$f(x_0)+\dfrac{x_n-x_0}{1!\,h}\,\Delta\,f(x_0)+\dfrac{(x_n-x_0)(x_n-x_1)}{2!\,h^2}\,\Delta^{2}\,f(x_0)+ $
  $+\dfrac{(x_n-x_0)(x_n-x_1)(x_n-x_2)}{3!\,h^3}\,\Delta^{3}\,f(x_0)+$
    $\overset{\underbrace{n+1}}{\ldots}+\dfrac{(x_n-x_0)(x_n-x_1)(x_n-x_2)\ldots(x_n-x_{n-1})}{n!\,h^n}\,\Delta^{n}\,f(x_0)$

En consecuencia, para un $x$ genérico, que no sea necesariamente una las abscisas de los datos, si bien su valor sea cercano a las mismas, será razonable escribir el polinomio interpolador de Newton de grado $n$ ( a partir de los $n+1$ puntos que se dan como datos ) como

$P_{n}\,(x)=f(x_0)+\dfrac{x-x_0}{1!\,h}\,\Delta\,f(x_0)+\dfrac{(x-x_0)(x-x_1)}{2!\,h^2}\,\Delta^{2}\,f(x_0)+ $
  $+\dfrac{(x-x_0)(x-x_1)(x-x_2)}{3!\,h^3}\,\Delta^{3}\,f(x_0)+$
    $\overset{\underbrace{n+1}}{\ldots}+\dfrac{(x-x_0)(x-x_1)(x-x_2)\ldots(x-x_{n-1})}{n!\,h^n}\,\Delta^{n}\,f(x_0)$

que, de manera abreviada, podemos escribir también de la forma $$\displaystyle P_{n}(x)=f(x_0)+\sum_{i=1}^{n}\,\dfrac{\Delta^{i}\,f(x_0)}{i!\,h^i}\,\prod_{j=0}^{i-1}\,(x-x_j)$$

-oOo-
Para el ejemplo concreto que nos hemos planteado ( recordemos la tabla de diferencias ):

podemos escribir el polinomio interpolador de Newton de grado $2$ ( para este conjunto de datos formado por tres puntos ):

$$P_{2}(x)=f(x_0)+\dfrac{x-x_0}{1! \cdot h}\,\Delta\,f(x_0)+\dfrac{(x-x_0)(x-x_1)}{2 !\cdot h^2}\,\Delta^{2}\,f(x_0)$$

con lo cual

$$P(x)=-1-\dfrac{x-(-2)}{1}\cdot (-1)+\dfrac{(x-(-2))(x-(-1))}{2}\cdot 2 = x^2+2x-1$$

Con ello, podríamos obtener la ordenada correspondiente a una abscisa próxima a alguna de las que figuran como datos; por ejemplo a un punto con abscisa $x=-3/2$, la ordenada que le corresponde, según el polinomio interpolador obtenido, es $P(-3/2)=(-/2)^2+2\cdot (3/2)-1=17/4$

$\square$

Obtención del polinomio interpolador por el método de Lagrange

En una artículo anterior había resuelto el siguiente problema de interpolación:

Sean tres puntos del plano cartesiano $A(-2,-1)$, $B(-1,-2)$ y $C(0,-1)$. Determínese el polinomio interpolador de segundo grado cuya gráfica pasa por estos tres puntos


Y lo había hecho calculando los coeficientes de la función polinómica interpoladora, que, al conocer $3$ puntos ha de ser de grado $3-1=2$: $f(x)=a_0+a_1\,x+a_2\,x^2$. El método que empleé, se basaba en resolver un sistema de ecuaciones lineales, lo cual es bastante laborioso en cuánto al número de operaciones aritméticas a realizar cuando lo generalizamos a $n+1$ puntos ( obteniendo una función polinómica interpoladora de grado $n$ ) para valores grandes de $n$, pues hay que resolver un sistema de de $n+1$ ecuaciones con $n+1$ incógnitas ( los coeficientes ). Voy a presentaros ahora otro método alternativo, el método de Lagrange, que es bastante más eficiente e ingenioso que el anterior.

Con los tres puntos como datos:

x     |   f(x)
-------------------
x_{0} |  f(x_0)
x_{1} |  f(x_1)
x_{2} |  f(x_2)

podemos escribir el polinomio interpolador de la forma $$f(x)=f(x_0)\,\Phi_{0}(x)+f(x_1)\,\Phi_{1}(x)+f(x_2)\,\Phi_{2}(x)$$ de manera que:

$$\Phi_{0}(x)=\left\{\begin{matrix}1 & \text{si} & x=x_0 \\ 0 & \text{si} & x=x_1 & \text{o bien} & x=x_2 \end{matrix}\right.$$

$$\Phi_{1}(x)=\left\{\begin{matrix}1 & \text{si} & x=x_1 \\ 0 & \text{si} & x=x_0 & \text{o bien} & x=x_2 \end{matrix}\right.$$

$$\Phi_{2}(x)=\left\{\begin{matrix}1 & \text{si} & x=x_2 \\ 0 & \text{si} & x=x_0 & \text{o bien} & x=x_1 \end{matrix}\right.$$

Además, dichas funciones han de ser polinómicas, de grado $2$, por lo que definiremos éstas de la siguiente manera:

$$\Phi_{0}(x)=\dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}$$

$$\Phi_{1}(x)=\dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}$$

$$\Phi_{2}(x)=\dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$$

con lo cual se tiene que el polinomio interpolador para los tres puntos dados es

$f(x)=f(x_0)\,\dfrac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+f(x_1)\,\dfrac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}+$
    $+f(x_2)\,\dfrac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$

Evidentemente, el método de interpolación de Lagrange se generaliza sin dificultad para $n+1$ puntos dados, obteniendo el polonomio interpolador ( que es único, independientemente del método que empleemos ) de grado $n$ $$\displaystyle f(x)=\sum_{i=0}^{n}\, f(x_i)\,\Phi_{i}(x)$$ donde $$\Phi_i=\prod_{j\neq i} \dfrac{x-x_j}{x_i-x_j}\; \; \text{para}\; \; i=0,1,\ldots,n$$

-oOo-
Teniendo en cuenta ahora los datos concretos del problema:

x  |  f(x)
----------
-2 |  -1
-1 |  -2
0  |  -1

Obtenemos el polinomio interpolador
$f(x)=(-1)\cdot \dfrac{(x-(-1))(x-0)}{(-2-0)(-2-(-1))}+(-2)\cdot \dfrac{(x-(-2))(x-0)}{(-1-(-2))(-1-0)}+$
    $+(-1)\,\dfrac{(x-(-2))(x-(-1))}{(0-(-2))(0-(-1))}=\ldots$
      $=x^2+2x-1$

Recordemos una vez más que el polinomio interpolador que pasa por un conjunto de $n+1$ puntos es de grado $n$ y único ( independientemente del método empleado para encontrarlo ); puede comprobarse que es el mismo que el obtenido por el método del cálculo directo de los coeficientes.

$\square$

Método básico para el cálculo del polinomio interpolador de grado n ( conocidas las coordenadas de un conjunto de n+1 puntos )

Consideremos que tenemos $n+1$ valores de una función polinómica de grado $n$, $f(x)=a_0+a_a\,x+a_2\,x^2+\ldots+a_n\,x^n$, cuyos coeficientes $a_0,a_1,\ldots,a_n$ son desconocidos. Y supongamos que conocemos $n+1$ puntos de la gráfica de la función de $f(x)$: $(x_i,f(x_i)$ para $i=0,1,2,\ldots,n$. Entonces podemos calcular los coeficientes imponiendo que la función satisfaga las siguientes $n+1$ condiciones $$\left\{ \begin{matrix}a_0+a_1\,x_0+a_2\,x_{0}^2+\ldots+a_{n+1}\,x_{0}^{n+1}&=&f(x_0) \\ a_0+a_1\,x_1+a_2\,x_{1}^2+\ldots+a_{n+1}\,x_{1}^{n+1}&=&f(x_1) \\ \ldots \\ a_0+a_1\,x_{n+1}+a_2\,x_{n+1}^2+\ldots+a_{n+1}\,x_{n+1}^{n+1}&=&f(x_{n+1}) \end{matrix}\right.$$ Así pues, para encontrar el polinomio interpolador ( que es único ), debemos resolver este sistema de $n+1$ ecuaciones lineales con $n+1$ incógnitas ( los coeficientes del sistema ).

Escribiendo el sistema en forma matricial $$A\,X=B$$ donde $A$ es la matriz $ (n+1) \times (n+1) $de los coeficientes del sistema es $$A=\begin{pmatrix}1 & x_0 & x_{0}^{2} & \ldots & x_{0}^{n}\\ 1 & x_1 & x_{1}^{2} & \ldots & x_{1}^{n} \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1 & x_n & x_{n}^{2} & \ldots & x_{n}^{n}\end{pmatrix}$$ $X$ es la matriz $n+1 \times 1$ de las incógnitas y $B$ es la matriz $n+1 \times 1$ de los términos independientes: $$X=\begin{pmatrix}a_0 \\ a_1 \\ a_2 \\ \ldots \\ a_n \end{pmatrix} \quad \quad \quad \quad B=\begin{pmatrix} f(x_0) \\ f(x_1) \\ f(x_2) \\ \ldots \\ f(x_n)\end{pmatrix}$$ Naturalmente, podemos resolverlo por Gauss, o por cualquier otro de los métodos estudiados; como, por ejemplo, por Cramer: $$a_0=\dfrac{\text{det}(A_0)}{\text{det}(A)}\;,\; a_1=\dfrac{\text{det}(A_1)}{\text{det}(A)}\;,\;a_2=\dfrac{\text{det}(A_2)}{\text{det}(A)}\;,\ldots,\;a_n=\dfrac{\text{det}(A_n)}{\text{det}(A)}$$ en el que reconocemos en su forma el determinante de Vandermonde $$\text{det}(A)=\begin{vmatrix}1 & a_0 & a_{0}^{2} & \ldots & a_{0}^{n}\\ 1 & a_1 & a_{1}^{2} & \ldots & a_{1}^{n} \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1 & a_n & a_{n}^{2} & \ldots & a_{n}^{n}\end{vmatrix}$$
siendo
$$\text{det}(A_0)=\begin{vmatrix}f(x_0) & a_0 & a_{0}^{2} & \ldots & a_{0}^{n}\\ f(x_1) & a_1 & a_{1}^{2} & \ldots & a_{1}^{n} \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ f(x_n) & a_n & a_{n}^{2} & \ldots & a_{n}^{n}\end{vmatrix}$$
$$\text{det}(A_1)=\begin{vmatrix}1 & f(x_0) & a_{0}^{2} & \ldots & a_{0}^{n}\\ 1 & f(x_1) & a_{1}^{2} & \ldots & a_{1}^{n} \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1 & f(x_2) & a_{n}^{2} & \ldots & a_{n}^{n}\end{vmatrix}$$
$$\text{det}(A_2)=\begin{vmatrix}1 & a_0 & f(x_0) & \ldots & a_{0}^{n}\\ 1 & a_1 & f(x_1) & \ldots & a_{1}^{n} \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1 & a_n & f(x_n) & \ldots & a_{n}^{n}\end{vmatrix}$$
$$\vdots$$
$$\text{det}(A_n)=\begin{vmatrix}1 & a_0 & a_{0}^{2} & \ldots & f(x_0)\\ 1 & a_1 & a_{1}^{2} & \ldots & f(x_1) \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 1 & a_n & a_{n}^{2} & \ldots & f(x_n)\end{vmatrix}$$

EJEMPLO. Sean tres puntos del plano cartesiano $A(-2,-1)$, $B(-1,-2)$ y $C(0,-1)$. Se pide:
a) Determinar el polinomio interpolador de segundo grado cuya gráfica pasa por estos tres puntos
b) Considerando un punto $P$ cuya abscisa es igual a $\dfrac{2}{3}$, ¿ cuál es la ordenada que le corresponde por interpolación ?

SOLUCIÓN.
a) La función interpoladora ( polinomio interpolador, de segundo grado en este caso, ya que contamos con 2+1 puntos de su gráfica ) tiene la forma $f(x)=a_0+a_1\,x+a_2\,x^2$. Para determinar el valor de los coeficientes $a_0,a_1$ y $a_2$, impondremos que la ecuación se cumpla para cada uno de los puntos ( método de los coeficientes indeterminados ). Así:
$$\left.\begin{matrix}A:&a_0+a_1\cdot (-2)+a_2 \cdot (-2)^2&=&-1 \\ B:&a_0+a_1\cdot (-1)+a_2 \cdot (-1)^2&=&-2 \\ C:&a_0+a_1\cdot 0+a_2 \cdot 0^2&=&-1 \end{matrix}\right\}$$ esto es
$$\left.\begin{matrix}a_0&+&-2\,a_1&+&4\,a_2 &=&-1 \\ a_0&-&a_1&+&a_2&=&-2 \\ a_0&&&&&=&-1 \end{matrix}\right\}$$

Entonces, empleando alguno de los métodos conocidos, pongamos que el de Cramer ( para ilustrar o dicho arriba, aunque en este caso tan sencillo, sería incluso más rápido si empleásemos el de Gauss ):

$$a_0=\dfrac{\begin{vmatrix}-1&-2&4\\-2&-1&1\\-1&0&0\end{vmatrix}}{\begin{vmatrix}1&-2&4\\1&-1&1\\1&0&0\end{vmatrix}}=-1$$

$$a_1=\dfrac{\begin{vmatrix}-1&-1&4\\1&-2&1\\1&-1&0\end{vmatrix}}{\begin{vmatrix}1&-2&4\\1&-1&1\\1&0&0\end{vmatrix}}=2$$

$$a_2=\dfrac{\begin{vmatrix}-1&-2&-1\\-2&-1&-2\\-1&0&-1\end{vmatrix}}{\begin{vmatrix}1&-2&4\\1&-1&1\\1&0&0\end{vmatrix}}=1$$

Por consiguiente, el polinomio interpolador es $$f(x)=x^2+2x-1$$

El siguiente gráfico muestra los tres puntos dados y la función que acabamos de deducir:

b)
Sustituyendo la variable independiente por el valor dado,
$$f(2/3)=\left(\dfrac{2}{3}\right)^2+2\cdot \dfrac{2}{3}-1$$
y operando llegamos al valor pedido $$f(2/3)=\dfrac{7}{9}$$ luego las coordenadas del punto pedido son $$P\left( \dfrac{2}{3},\dfrac{7}{9}\right)$$
$\square$

Ejemplo de interpolación por el método de Lagrange haciendo uso de MAXIMA

[nota del autor]