Graficación

Antes que nada, tenemos que aprender a graficar en Python, lo manera mas fácil de graficar es usando la función plot de la libería matplotlib, asi que importamos esta función:


In [ ]:
from matplotlib.pyplot import plot

y la usamos como cualquier función, dandole dos listas, una con todos los valores de $x$, y otra con todos los valores de $y$:


In [ ]:
plot([0,1], [2,3])

Sin embargo no aparece nada, para esto es necesario decirle a la librería matplotlib que queremos que nos muestre las gráfcias en linea con nuestro notebook, por lo que utilizamos:


In [ ]:
%matplotlib inline

y si graficamos ahora:


In [ ]:
plot([0,1],[2,3])

De la misma manera, podemos gráficar cualquier función, creada o importada, empecemos con:

$$ y=\sin{x} $$

Antes que nada, importemos esta función de la librería numpy, así como la función linspace:


In [ ]:
from numpy import sin

In [ ]:
from numpy import linspace

La función linspace nos ayudara a crear un arreglo lineal de datos el cual define al eje $x$:


In [ ]:
xs = linspace(0, 10, 100)

In [ ]:
xs

Ahora solo tenemos que meter estos datos a la función $\sin$:


In [ ]:
ys = sin(xs)

In [ ]:
ys

Con lo que vemos que tenemos nuestros datos, pero lo principal es graficarlo, por lo que le damos estos dos arreglos a plot y obtendremos:


In [ ]:
plot(xs, ys)

Una opción que puedo utilizar para darle formato a mi gráfica es "o", con lo que nos mostrará los datos como puntos, en lugar de una linea que conecta a estos:


In [ ]:
plot(xs, ys, "o")

y despues de esta breve introducción a la graficación en Python, empecemos con nuestro problema:

Lo que queremos es una función que pase exactamente por los siguientes puntos:

$i$ 0 1 2 3
$f(x_i)$ -3 0 5 7
$x_i$ 0 1 3 6

El primer paso será guardar estos datos en variables de Python, especificamente listas:


In [ ]:
datos_x = [0, 1, 3, 6]
datos_y = [-3, 0, 5, 7]

Utilizando la formula de interpolación por polinomios de Lagrange, tenemos que:

$$ p(x) = L_0(x)f(x_0) + L_1(x)f(x_1) + L_2(x)f(x_2) $$

en donde cada uno de los polinomios de Lagrange, se calcula con la formula:

$$ L_i(x) = \prod_{j=0, j\ne i}^n \frac{x-x_j}{x_i-x_j} $$

en una sola linea, esto se ve asi:


In [ ]:
L0 = lambda x: ((x - datos_x[1])*(x - datos_x[2])*(x - datos_x[3]))/((datos_x[0] - datos_x[1])*(datos_x[0] - datos_x[2])*(datos_x[0] - datos_x[3]))

In [ ]:
L0(5)

Sin embargo esta solución solo calcula uno de los polinomios de Lagrange y solo sirve para el caso de $4$ datos; el primer paso para crear una solución total, es separar estos polinomios en cada uno de los multiplicandos:


In [ ]:
L01 = lambda x: (x - datos_x[1])/(datos_x[0] - datos_x[1])
L02 = lambda x: (x - datos_x[2])/(datos_x[0] - datos_x[2])
L03 = lambda x: (x - datos_x[3])/(datos_x[0] - datos_x[3])

L0  = lambda x: L01(x)*L02(x)*L03(x)

In [ ]:
L01(5)

In [ ]:
L02(5)

In [ ]:
L03(5)

In [ ]:
L0(5)

Por lo que vemos que este método es equivalente, tan solo tenemos que crear un metodo iterativo para encontrar estos polinomios.

Quiero hacer notar que estos objetos que creamos son funciones, de la misma manera que sus elementos:


In [ ]:
L0

In [ ]:
L01, L02, L03

Tambien notar, que para la creación de estas funciones, utilice un subconjunto de estos datos, que excluia al primer elemento; crearé una lista que excluya justo a este primer elemento:


In [ ]:
dato_x0 = datos_x[0]
datos_x0 = datos_x[1:]

In [ ]:
dato_x0

In [ ]:
datos_x0

y con un ciclo for, agregaré funciones en una lista, en la cual estarán cada uno de los multiplicandos del primer polinomio de Lagrange:


In [ ]:
L0s = []
for i in range(len(datos_x0)):
    L0s.append(lambda x, i=i: (x - datos_x0[i])/(dato_x0 - datos_x0[i]))

Note que L0s esta compuesto de funciones, exactamente igual que L01, L02 y L03:


In [ ]:
L0s

Ahora solo tengo que multiplicar estas funciones, para obtener L0, utilizaré una función llamada reduce, la cual evaluará por pares, segun la regla que le doy (multiplicación):


In [ ]:
from functools import reduce

In [ ]:
L0 = reduce(lambda x, y: lambda z :x(z)*y(z), L0s)

In [ ]:
L0(5)

Una vez que tengo L0, puedo obtener los siguientes polinomios:


In [ ]:
dato_x1 = datos_x[1]
datos_x1 = datos_x[:1] + datos_x[2:]

In [ ]:
L1s = []
for i in range(len(datos_x1)):
    L1s.append(lambda x, i=i: (x - datos_x1[i])/(dato_x1 - datos_x1[i]))

In [ ]:
L1 = reduce(lambda x, y: lambda z :x(z)*y(z), L1s)

In [ ]:
L1(5)

Pero esto me obligaría a escribir código tantas veces como datos tenga, lo cual no es aceptable; de la misma manera en que los multiplicandos los metí en una lista, agregaré estos polinomios de Lagrange en una lista:


In [ ]:
Ls = []
for j in range(len(datos_x)):
    dato_xi = datos_x[j]
    datos_xi = datos_x[:j] + datos_x[j+1:]
    
    Lis = []
    for i in range(len(datos_xi)):
        Lis.append(lambda x, i=i, dato_xi=dato_xi, datos_xi=datos_xi: (x - datos_xi[i])/(dato_xi - datos_xi[i]))
        
    Li = reduce(lambda x, y: lambda z: x(z)*y(z), Lis)
    Ls.append(Li)

In [ ]:
Ls

Con lo que tenemos las funciones asociadas a cada uno de los polinomios de Lagrange, estos a su vez, tienen que se multiplicados por los datos_y, por lo que creamos para el primer polinomio:


In [ ]:
Lf0 = lambda x: Ls[0](x)*datos_y[0]
Lf0(5)

O bien, para todos de un solo jalón:


In [ ]:
Lfs = []
for j in range(len(datos_y)):
    Lfi = lambda x, j=j: Ls[j](x)*datos_y[j]
    Lfs.append(Lfi)

Por ultimo, todos estos terminos que se encuentran guardados en Lfs, podemos sumarlos usando reduce nuevamente, pero utilizando la regla de suma ahora:

$$ p(x) = \sum_{i=0}^n L_i(x) f(x_i) $$

In [ ]:
interp = reduce(lambda x, y: lambda z: x(z)+y(z), Lfs)

Y al evaluarlo en un número real, se comporta como esperamos:


In [ ]:
interp(5)

Si ahora creamos un arreglo de datos desde $0$ hasta $10$, tan solo para asegurar que vemos todos los datos:


In [ ]:
xs = linspace(0, 10, 100)
ys = interp(xs)

In [ ]:
plot(xs, ys)

Y graficar esta función en conjunto con los datos originales es facil:


In [ ]:
plot(xs, ys)
plot(datos_x, datos_y, "o")

De tal manera que podemos crear una función que haga todo el trabajo por nosotros:


In [ ]:
def interpolacion_Lagrange(datos_x, datos_y):
    Ls = []
    
    for j in range(len(datos_x)):
        dato_xi = datos_x[j]
        datos_xi = datos_x[:j] + datos_x[j+1:]

        Lis = []
        for i in range(len(datos_xi)):
            Lis.append(lambda x, i=i, dato_xi=dato_xi, datos_xi=datos_xi: (x - datos_xi[i])/(dato_xi - datos_xi[i]))

        Li = reduce(lambda x, y: lambda z: x(z)*y(z), Lis)
        Ls.append(Li)
        
    Lfs = []
    for j in range(len(datos_y)):
        Lfi = lambda x, j=j: Ls[j](x)*datos_y[j]
        Lfs.append(Lfi)
        
    interp = reduce(lambda x, y: lambda z: x(z)+y(z), Lfs)
    
    return interp

In [ ]:
dx = [0, 1, 3, 6]
dy = [-3, 0, 5, 7]
poli = interpolacion_Lagrange(dx, dy)

In [ ]:
xs = linspace(0, 6, 100)
ys = poli(xs)

In [ ]:
plot(xs, ys)
plot(dx, dy , "o")

Problemas

  1. Obtenga $6$ datos al azar, en el espacio $\Re^{2}$ y obtenga la interpolación polinómica por polinomios de Largange.
  2. Resuelva el problema 4.6 del libro de texto, ubicado en la pagina 387