Interpolación polinómica

Contenidos

Polinomio interpolante

Tenemos una serie de puntos $(x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)$ y queremos encontrar un polinomio que pase por todos ellos.

Sea el conjunto de puntos: $$C=\{(2,2),(3,6),(4,5),(5,5),(6,6)\}$$

Vamos a dibujar los puntos.

Para que inserte las figuras en la notebook, en lugar de abrir cada figura en una ventana nueva.

In [2]:
%matplotlib inline

Importamos los módulos matplotlib.pyplot y numpy.

In [3]:
import numpy as np
import matplotlib.pyplot as plt

Los puntos los almacenaremos en matrices

In [4]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Y los dibujamos

In [5]:
plt.plot(x,y,'ro')
plt.axis([1.5,6.5,1.5,6.5]);   # Dibujamos en [1.5,6.5]x[1.5,6.5]

Solución con la matriz de Vandermonde

Como tenemos $5$ puntos podemos plantear $5$ ecuaciones y necesitamos $5$ incógnitas que serán los coeficientes del polinomio. Por lo tanto:

$$P_4(x)=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4$$

es decir $a_0,a_1,a_2,a_3,a_4$, que son $5$ incógnitas y el polinomio es, en principio, de grado $4$. En general, para los puntos $(x_0,y_0),(x_1,y_1),\ldots,(x_n,y_n)$ el polinomio interpolante será de grado $\leq n$.

Las ecuaciones serían:

$$ \begin{array}{lll} P_{4}(x_{0})=y_{0} & \quad & a_{0}+a_{1}x_{0}+a_{2}x_{0}^{2}+a_{3}x_{0}^{3}+a_{4}x_{0}^{4}=y_{0}\\ P_{4}(x_{1})=y_{1} & \quad & a_{0}+a_{1}x_{1}+a_{2}x_{1}^{2}+a_{3}x_{1}^{3}+a_{4}x_{1}^{4}=y_{1}\\ P_{4}(x_{2})=y_{2} & \quad & a_{0}+a_{1}x_{2}+a_{2}x_{2}^{2}+a_{3}x_{2}^{3}+a_{4}x_{2}^{4}=y_{2}\\ P_{4}(x_{3})=y_{3} & \quad & a_{0}+a_{1}x_{3}+a_{2}x_{3}^{2}+a_{3}x_{3}^{3}+a_{4}x_{3}^{4}=y_{3}\\ P_{4}(x_{4})=y_{4} & \quad & a_{0}+a_{1}x_{4}+a_{2}x_{4}^{2}+a_{3}x_{4}^{3}+a_{4}x_{4}^{4}=y_{4} \end{array} $$

Y el sistema a resolver es:

$$ \left(\begin{array}{ccccc} 1 & x_{0} & x_{0}^{2} & x_{0}^{3} & x_{0}^{4}\\ 1 & x_{1} & x_{1}^{2} & x_{1}^{3} & x_{1}^{4}\\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3} & x_{2}^{4}\\ 1 & x_{3} & x_{3}^{2} & x_{3}^{3} & x_{3}^{4}\\ 1 & x_{4} & x_{4}^{2} & x_{4}^{3} & x_{4}^{4} \end{array}\right)\left(\begin{array}{c} a_{0}\\ a_{1}\\ a_{2}\\ a_{3}\\ a_{4} \end{array}\right)=\left(\begin{array}{c} y_{0}\\ y_{1}\\ y_{2}\\ y_{3}\\ y_{4} \end{array}\right) $$

La matriz de coeficientes del sistema se llama Matriz de Vandermonde.

Ejercicio 1

  • Escribir una función Vandermonde(x) que tenga como argumentos de entrada el vector x y como argumento de salida la matriz de Vandermonde V.
  • Calcular los coeficientes del polinomio resolviendo el sistema Va=y. Usar la orden numpy.linalg.solve.
  • Dibujar el polinomio interpolante con los nodos. Utilizar la orden numpy.polyval para evaluar el polinomio en un punto o puntos a partir de los coeficientes obtenidos en el apartado anterior.

Usar los nodos:

In [6]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Comprobar que el programa también funciona con los nodos:

In [7]:
x = np.array([0.,1.,2.,3.,4.,5.,6.])
y = np.array([3.,5.,6.,5.,4.,4.,5.])
In [8]:
%run Ejercicio1.py

Como resultados intermedios en este problema, podemos obtener la matriz de Vandermonde y los coeficientes del polinomio.

In [9]:
# opciones del print
np.set_printoptions(precision = 2)   # solo dos decimales
np.set_printoptions(suppress = True) # no usar notación exponencial

print 'Matriz de Vandermonde V'
print V

Matriz de Vandermonde V
[[    1.     2.     4.     8.    16.]
 [    1.     3.     9.    27.    81.]
 [    1.     4.    16.    64.   256.]
 [    1.     5.    25.   125.   625.]
 [    1.     6.    36.   216.  1296.]]
In [10]:
print 'Coeficientes del polinomio (desde a0 hasta a4)'
print a

Coeficientes del polinomio (desde a0 hasta a4)
[-75.    81.   -29.25   4.5   -0.25]

Existe una función Python que calcula la Matriz de Vandermonde. El orden de los elementos es distinto al que hemos usado

In [11]:
va = np.vander(x, len(x))
print(va)

[[   16.     8.     4.     2.     1.]
 [   81.    27.     9.     3.     1.]
 [  256.    64.    16.     4.     1.]
 [  625.   125.    25.     5.     1.]
 [ 1296.   216.    36.     6.     1.]]

Pero este no es un método aconsejable porque la matriz de Vandermonde, en general, está mal condicionada y puede dar lugar a grandes errores en la solución del sistema.

Interpolación mediante los polinomios fundamentales de Lagrange

Para cada $k=0,1,\ldots,n$, existe un único polinomio $\ell_{k}$ de grado $\leq n$ tal que $\ell_{k}\left(x_{j}\right)=\delta_{kj}$ (uno si $k=j$, cero en caso contrario):

$$ \ell_{k}\left(z\right)= \frac{(z-x_{0})\ldots(z-x_{k-1})(z-x_{k+1})\ldots(z-x_{n})}{(x_{k}-x_{0}) \ldots(x_{k}-x_{k-1})(x_{k}-x_{k+1})\ldots(x_{k}-x_{n})}. $$

Los polinomios $\ell_{0},\ell_{1},\ldots,\ell_{n}$ son los polinomios fundamentales de Lagrange.

El polinomio interpolante de Lagrange en los puntos $x_{0,}x_{1},\ldots,x_{n}$ relativo a los valores $y_{0},y_{1},\ldots,y_{n}$ es

$$ P_{n}\left(x\right)=y_{0}\ell_{0}\left(x\right)+ y_{1}\ell_{1}\left(x\right)+\cdots+y_{n}\ell_{n}\left(x\right) $$

Ejercicio 2

Escribir un programa que calcule el polinomio interpolante de Lagrange utilizando los polinomios funcamentales de Lagrange. Para ello, crear dos funciones. La primera será lagrange_fund(k,z,x) que tendrá:

  • Argumentos de entrada: índice $k$ del polinomio Fundamental, un punto (o vector) de evaluación, $z$ y los nodos $x$ del problema de interpolación.
  • Argumentos de salida: el valor del $k-$ésimo polinomio fundamental correspondiente a los nodos en el punto (o vector) $z$.

y una funcion polinomio_lagrange(z,x,y) que tendrá:

  • Argumentos de entrada: un punto (o vector) a evaluar $z$ y los nodos $x$, $y$ del problema.
  • Argumentos de salida: el valor del polinomio interpolante de Lagrange en $z$.

Dibujar el polinomio Fundamental de Lagrange con índice $k=2$ en el intervalo $[2,6]$. Dibujar el polinomio interpolante en el intervalo $[2,6]$.

Usar los nodos:

In [12]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Comprobar que también funciona con los nodos:

In [13]:
x = np.array([0.,1.,2.,3.,4.,5.,6.])
y = np.array([3.,5.,6.,5.,4.,4.,5.])
In [14]:
%run Ejercicio2.py

El inconveniente de este método es que si queremos incorporar un nodo nuevo tenemos que rehacer todos los cálculos.

Interpolación mediante diferencias divididas

Fórmula de Newton

El polinomio de interpolación de Lagrange en los puntos $x_{0,}x_{1},\ldots,x_{n}$ relativo a los valores $y_{0},y_{1},\ldots,y_{n}$ es

$$ P_{n}\left(x\right)= \left[y_{0}\right]+\left[y_{0},y_{1}\right]\left(x-x_{0}\right)+ \left[y_{0},y_{1},y_{2}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots +\left[y_{0},y_{1},\ldots,y_{n}\right]\left(x-x_{0}\right)\left(x-x_{1}\right)\cdots\left(x-x_{n-1}\right) $$

Tabla de diferencias divididas

Los coeficientes del polinomio anterior son la primera fila de la tabla: $$ \begin{array}{cccccc} x_{0} & y_{0} & \left[y_{0},y_{1}\right]=\frac{y_{0}-y_{1}}{x_{0}-x_{1}} & \left[y_{0},y_{1},y_{2}\right]=\frac{\left[y_{0},y_{1}\right]-\left[y_{1},y_{2}\right]}{x_{0}-x_{2}} & \ldots & \left[y_{0},y_{1},\ldots,y_{n}\right]=\frac{\left[y_{0},\ldots,y_{n-1}\right]-\left[y_{1},\ldots,y_{n}\right]}{x_{0}-x_{n}}\\ x_{1} & y_{1} & \left[y_{1},y_{2}\right]=\frac{y_{1}-y_{2}}{x_{1}-x_{2}} & \left[y_{1},y_{2},y_{3}\right]=\frac{\left[y_{1},y_{2}\right]-\left[y_{2},y_{3}\right]}{x_{1}-x_{3}} & \ldots\\ x_{2} & y_{2} & \left[y_{2},y_{3}\right]=\frac{y_{2}-y_{3}}{x_{2}-x_{3}} & \left[y_{2},y_{3},y_{4}\right]=\frac{\left[y_{2},y_{3}\right]-\left[y_{3},y_{4}\right]}{x_{2}-x_{4}} & \ldots\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ x_{n-1} & y_{n-1} & \left[y_{n-1},y_{n}\right]=\frac{y_{n-1}-y_{n}}{x_{n-1}-x_{n}}\\ x_{n} & y_{n} \end{array} $$

Con este método, si añadimos un nodo, no hay que rehacer todos los cálculos, sino que se añade una línea más a la tabla anterior en su parte inferior.

Ejercicio 3

Escribir un programa que calcule el polinomio interpolante de Lagrange utilizando las diferencias divididas. Para ello, crear dos funciones. La primera será dif_div(x,y) que tendrá:

  • Argumentos de entrada: los nodos $x$, $y$ del problema de interpolación.
  • Argumentos de salida: una matriz que contenga las diferencias divididas

y una funcion polinomio_newton(z,x,y) que tendrá:

  • Argumentos de entrada: un punto (o vector) a evaluar $z$ y los nodos $x$, $y$ del problema.
  • Argumentos de salida: el valor del polinomio interpolante de Lagrange en $z$.

Escribir la matriz de diferencias divididas. Dibujar el polinomio interpolante en el intervalo $[2,6]$.

Usar los nodos:

In [15]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])

Comprobar que también funciona con los nodos:

In [16]:
x = np.array([0.,1.,2.,3.,4.,5.,6.])
y = np.array([3.,5.,6.,5.,4.,4.,5.])
In [17]:
%run Ejercicio3.py

[[ 2.    2.    4.   -2.5   1.   -0.25]
 [ 3.    6.   -1.    0.5   0.    0.  ]
 [ 4.    5.    0.    0.5   0.    0.  ]
 [ 5.    5.    1.    0.    0.    0.  ]
 [ 6.    6.    0.    0.    0.    0.  ]]

Interpolación con funciones de python

El polinomio de interpolación se puede obtener con la función de numpy polyfit escogiendo como grado del polinomio el número de puntos menos uno. Nos da los coeficientes del polinomio interpolante en un vector, empezando por el de mayor grado:

In [18]:
x = np.array([2.,3.,4.,5.,6.])
y = np.array([2.,6.,5.,5.,6.])
pol = np.polyfit(x,y,len(x)-1)  # coeficientes del polinomio

xx = np.linspace(min(x),max(x))
yy = np.polyval(pol,xx)         # valor del polinomio en los puntos de la matriz xx

Representamos el polinomio

In [19]:
plt.plot(xx, yy, '-',x, y, 'ro')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Ejercicio 4

Escribir una función chebyshev(f,a,b,n) que interpole la función f en el intervalo [a,b] utilizando n nodos:

  • Utilizando nodos equiespaciados que incluyan los extremos del intervalo.
  • Utilizando nodos de Chebyshev, que se definen en el intervalo $[-1,1]$ según la fórmula:
$$ x_{i}^{(n)}=\cos\frac{\left(2i-1\right)\pi}{2n}\quad i=1,2,\ldots,n $$

El número de nodos es $n = 11$, el itervalo $[a,b]=[-1,1]$ y la función $f$

$$ f(x)=\frac{1}{1+25x^{2}} $$

Para interpolar utilizar las funciones numpy.polyfit y numpy.polyval.

En ambos casos dibujar los nodos, la función y el polinomio interpolador.

Comprobar que también funciona con 15 nodos y con la función:

$$ f(x)=e^{-20x^2} $$
In [20]:
%run Ejercicio4.py

Podemos obtener la interpolación a trozos con funciones de scipy.

In [21]:
from scipy.interpolate import interp1d

f1 = interp1d(x, y, kind='nearest')
f2 = interp1d(x, y, kind='linear')
f3 = interp1d(x, y, kind='cubic')

xx = np.linspace(min(x),max(x))

De grado cero

In [22]:
plt.plot(x, y, 'ro', xx, f1(xx), '-')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

De grado uno o lineal

In [23]:
plt.plot(x, y, 'ro', xx, f2(xx), '-')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

De grado tres o cúbica

In [24]:
plt.plot(x, y, 'ro', xx, f3(xx), '-')
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Y la interpolación con splines (con la opción not-a-knot) con una función de scipy

In [25]:
from scipy import interpolate
In [26]:
s = interpolate.InterpolatedUnivariateSpline(x, y)
xx = np.linspace(min(x),max(x))
yy = s(xx)

plt.figure()
plt.plot(x, y, 'ro', xx, yy)
plt.axis([min(xx)-d, max(xx)+d, min(yy)-d, max(yy)+d]);

Ejercicio 5

Realizar una programa que, para la función $f(x) = \cos x$ en el intervalo $[0,10]$ tomando como nodos $x=0,1,\ldots,10$ dibuje los polinomios de interpolación:

  • Lineal a trozos.
  • Con splines cúbicas.

y calcule, para los puntos

xx = np.linspace(0,10,100)

el error absoluto en cada caso. Para evaluar el error utilizar la función numpy.linalg.norm.

$$E_a = \left\Vert P(xx)-f(xx) \right\Vert $$
In [27]:
%run Ejercicio5.py

Error interpolación lineal a trozos = 0.64506
Error interpolación con splines     = 0.06776