Página web del curso

Interpolación polinómica

Importamos los módulos matplotlib.pyplot, numpy y las funciones de numpy para polinomios.

In [22]:
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as pol

Polinomio interpolante de Lagrange

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=\{(-1,1),(0,3),(2,4),(3,3),(5,1)\}$$

Los puntos los almacenaremos en matrices numpy (numpy arrays)

In [13]:
x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Y dibujamos los puntos y el polinomio de Lagrange, que es el polinomio que pasa por todos ellos

In [14]:
xp = np.linspace(min(x),max(x))
p  = pol.polyfit(x,y,len(x)-1)
yp = pol.polyval(xp,p)

plt.figure()
plt.plot(xp,yp,'b-', label = 'polinomio interpolante')
plt.plot( x, y,'ro', label = 'puntos')
plt.legend()
plt.show()

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})}{(x_{k}-x_{0})}\cdots \frac{(z-x_{k-1})}{(x_{k}-x_{k-1})} \cdot \frac{(z-x_{k+1})}{(x_{k}-x_{k+1})} \cdots \frac{(z-x_{n})}{(x_{k}-x_{n})} $$

Por ejemplo, con nuestros 5 nodos, serían

$$ \ell_{0}\left(z\right)= \Bigg\rvert^{x_0} \frac{(z-x_1)}{(x_0-x_1)} \cdot \frac{(z-x_2)}{(x_0-x_2)} \cdot \frac{(z-x_3)}{(x_0-x_3)} \cdot \frac{(z-x_4)}{(x_0-x_4)} $$$$ \ell_{1}\left(z\right)= \frac{(z-x_0)}{(x_1-x_0)} \Bigg\rvert^{x_1} \frac{(z-x_2)}{(x_1-x_2)} \cdot \frac{(z-x_3)}{(x_1-x_3)} \cdot \frac{(z-x_4)}{(x_1-x_4)} $$$$ \ell_{2}\left(z\right)= \frac{(z-x_0)}{(x_2-x_0)} \cdot \frac{(z-x_1)}{(x_2-x_1)} \Bigg\rvert^{x_2} \frac{(z-x_3)}{(x_2-x_3)} \cdot \frac{(z-x_4)}{(x_2-x_4)} $$$$ \ell_{3}\left(z\right)= \frac{(z-x_0)}{(x_3-x_0)} \cdot \frac{(z-x_1)}{(x_3-x_1)} \cdot \frac{(z-x_2)}{(x_3-x_2)} \Bigg\rvert^{x_3} \frac{(z-x_4)}{(x_3-x_4)} $$$$ \ell_{4}\left(z\right)= \frac{(z-x_0)}{(x_4-x_0)} \cdot \frac{(z-x_1)}{(x_4-x_1)} \cdot \frac{(z-x_2)}{(x_4-x_2)} \cdot \frac{(z-x_3)}{(x_4-x_3)} \Bigg\rvert^{x_4} $$

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(z\right)=y_{0}\ell_{0}\left(z\right)+ y_{1}\ell_{1}\left(z\right)+\cdots+y_{n}\ell_{n}\left(z\right) $$

Ejercicio 1

Escribir un programa que calcule los polinomios fundamentales de Lagrange. Para ello, crear una función lagrange_fund(k,x,z) que tendrá:

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

Dibujar todos los polinomios Fundamental de Lagrange.

Usar los nodos:

x = np.array([-1., 0, 2, 3, 5])

Nota:

  • Puedes empezar haciendo la función para un punto $z$ y luego modificarla para que admita un vector. Por ejemplo, la primera versión funcionaría
In [15]:
k = 2
z = 1.3
yz = lagrange_fundamental(k,x,z)
print(yz)

1.0448388888888889

y ahora sabemos que $\ell_{2}\left(1.3\right)=1.0448388888888889$. La segunda versión haría

In [16]:
k = 2
z = np.array([1.3, 2.1, 3.2])
yz = lagrange_fundamental(k,x,z)
print(yz)

[ 1.04483889  0.94395    -0.2688    ]

y $\ell_{2}\left([1.3, \,2.1,\,3.2]\right)=[ 1.04483889,\,0.94395,\, -0.2688 ]$.

  • Para dibujar los puntos, se puede tener en cuenta que sus $x$s están contenidas en el vector de nodos x y sus correspondientes $y$s se pueden obtener como las filas de la matriz identidad:
In [17]:
print(np.eye(len(x)))

[[1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0.]
 [0. 0. 1. 0. 0.]
 [0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 1.]]
In [18]:
%run Ejercicio1.py

Recordamos que 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(z\right)=y_{0}\ell_{0}\left(z\right)+ y_{1}\ell_{1}\left(z\right)+\cdots+y_{n}\ell_{n}\left(z\right) $$

Ejercicio 2

Escribir un programa que calcule el polinomio interpolante de Lagrange utilizando los polinomios fundamentales de Lagrange. Para ello, usar la función lagrange_fund(k,x,z) y una funcion polinomio_lagrange(x,y,z) que tendrá:

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

Dibujar el polinomio interpolante en el intervalo $[-1,5]$.

Usar los nodos:

x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Comprobar que también funciona con los nodos en el intervalo $[-1,7]$:

x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])
In [19]:
%run Ejercicio2.py

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

 Nodos de Chebyshev

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 [23]:
x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])
p = pol.polyfit(x,y,len(x)-1)  # coeficientes del polinomio

xp = np.linspace(min(x),max(x))
yp = pol.polyval(xp,p)         # valor del polinomio en los puntos de xp

Representamos el polinomio

In [24]:
plt.figure()
plt.plot(xp, yp, '-',x, y, 'o')
plt.show()

Ejercicio 3

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}=\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_1(x)=\frac{1}{1+25x^{2}} $$

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

$$ f_2(x)=e^{-20x^2} $$

Seguir los pasos:

  • Construir los nodos: los n nodos se almacenarán en un vector x de n elementos en ambos casos. Para los nodos equiespaciados se pueden usar np.linspace o np.arange. Las y correspondientes se obtienen utilizando f1 y f2.

  • Construir los polinomios interpolantes: Usar las funciones pol.polyfit y pol.polyval para calcular los polinomios interpolantes.

  • Dibujar la gráfica: En ambos casos dibujar los puntos, la función y el polinomio interpolador. Usar plt.axis([-1.05, 1.05, -0.3, 2.3]) para que dibuje todas las gráficas en el rectángulo $[-1.05, 1.05]\times[ -0.3, 2.3]$. Usar más de 50 puntos (valor por defecto) en np.linspace porque el polinomio interpolante oscila mucho.

In [25]:
%run Ejercicio3.py

------------  Función f1  ------------


------------  Función f2  ------------

Ejercicios propuestos

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(z\right)= \left[y_{0}\right]+\left[y_{0},y_{1}\right]\left(z-x_{0}\right)+ \left[y_{0},y_{1},y_{2}\right]\left(z-x_{0}\right)\left(z-x_{1}\right)+\cdots +\left[y_{0},y_{1},\ldots,y_{n}\right]\left(z-x_{0}\right)\left(z-x_{1}\right)\cdots\left(z-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} $$

Ejercicio 4

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(x,y,z) que tendrá:

  • Argumentos de entrada: los nodos x, y del problema y un punto (o vector) a evaluar z.
  • 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 de interpolación.

Usar los nodos:

x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Comprobar que también funciona con los nodos:

x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])
In [26]:
%run Ejercicio4.py

TABLA DE DIFERENCIAS DIVIDIDAS
[[-1.    1.    2.   -0.5   0.    0.02]
 [ 0.    3.    0.5  -0.5   0.1   0.  ]
 [ 2.    4.   -1.    0.    0.    0.  ]
 [ 3.    3.   -1.    0.    0.    0.  ]
 [ 5.    1.    0.    0.    0.    0.  ]]

Interpolación polinomial a trozos con funciones de python

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

In [27]:
from scipy.interpolate import interp1d, CubicSpline

f1 = interp1d(x, y, kind='nearest')
f2 = interp1d(x, y, kind='linear')
f3 = CubicSpline(x, y, bc_type='natural')

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

De grado cero

In [28]:
plt.figure()
plt.plot(x, y, 'o', xp, f1(xp), '-')
plt.show()

De grado uno o lineal

In [29]:
plt.figure()
plt.plot(x, y, 'o', xp, f2(xp), '-')
plt.show()

De grado tres, un spline cúbico natural

In [30]:
plt.figure()
plt.plot(x, y, 'o', xp, f3(xp), '-')
plt.show()

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

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

el error absoluto en cada caso. Para evaluar el error utilizar la función Ea = np.linalg.norm(yp-f(xp)).

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

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

Solución con la matriz de Vandermonde

Para los puntos

In [33]:
x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

si queremos calcular el polinomio que pasa por todos ellos, 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) $$

es decir, $Vp=y$, donde $V$ es la matriz de coeficientes del sistema, Matriz de Vandermonde, $y$ contiene los términos independientes del sistema y la solución del sistema son los coeficientes del polinomio, contenidos en $p$.

Ejercicio 6

  1. Escribir una función Vandermonde(x) que tenga como argumento de entrada el vector x y como argumento de salida la matriz de Vandermonde V.
  2. Escribir una función polVandermonde(x,y) que tenga como argumento de entrada los valores de las coordenadas x e y de los nodos y como salida p, que contiene los coeficientes del polinomio interpolante que pasa por los nodos, es decir, $p = (a_0,a_1,a_2,a_3,a_4)$. Para ello:
    • Calcular V llamando a la función Vandermonde.
    • Calcular los coeficientes del polinomio resolviendo el sistema lineal Vp = y. Usar la orden p = np.linalg.solve(V,y).
  3. Dibujar el polinomio interpolante con los nodos:
    • xp = np.linspace(min(x),max(x)) nos da 50 puntos $x$ en el intervalo de interpolación.
    • Utilizar la orden yp = pol.polyval(xp,p) para evaluar el polinomio en un punto o puntos xp a partir de los coeficientes del polinomio, p, obtenidos en el apartado anterior y obtener sus valores yp.

Usar los nodos:

x = np.array([-1., 0, 2, 3, 5])
y = np.array([ 1., 3, 4, 3, 1])

Comprobar que también funciona con los nodos:

x1 = np.array([-1., 0, 2, 3, 5, 6, 7])
y1 = np.array([ 1., 3, 4, 3, 2, 2, 1])
In [32]:
%run Ejercicio6.py

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

In [33]:
print('Matriz de Vandermonde V')
print(V)

Matriz de Vandermonde V
[[  1.  -1.   1.  -1.   1.]
 [  1.   0.   0.   0.   0.]
 [  1.   2.   4.   8.  16.]
 [  1.   3.   9.  27.  81.]
 [  1.   5.  25. 125. 625.]]
In [34]:
print('Coeficientes del polinomio')
print(p)

Coeficientes del polinomio
[ 3.    1.6  -0.48 -0.07  0.02]

Existe una función numpy que calcula la Matriz de Vandermonde

In [37]:
va = pol.polyvander(x, len(x)-1)
print(va)

[[  1.  -1.   1.  -1.   1.]
 [  1.   0.   0.   0.   0.]
 [  1.   2.   4.   8.  16.]
 [  1.   3.   9.  27.  81.]
 [  1.   5.  25. 125. 625.]]

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.