GT_MM_1718_T3_Tutorial

587 days ago by GT_MM_1718

Interpolación polinomial

Método de Lagrange

Vamos a calcular el polinomio interpolador usando el método de Lagrange.

En primer lugar construiremos una función que calcula los polinomios de Lagrange:

def lagrange(px,i): n=len(px)-1 return prod([ x-px[k] for k in [0..n] if k!=i ])/prod([px[i]-px[k] for k in [0..n] if k!=i]) 
       

La probamos calculando el primer polinomio de Lagrange de los puntos -1,1,3.

plot(lagrange([-1,1,3],0),(x,-1,3)) 
       

Vamos el polinomio interpolador de (-1,3), (1,2), (3,-1).

Para calcular el polinomio interpolador basta multiplicar cada polinomio de Lagrange por la coordenada $y$ correspondiente y sumar el resultado:

px=[-1,1,3] p(x)=3*lagrange(px,0)+2*lagrange(px,1)+-1*lagrange(px,2) points([(-1,3),(1,2),(3,-1)])+plot(p(x),(x,-1,3)) 
       

Diferencias divididas o método de interpolación de Newton

La función diferencias_divididas, a continuación, devuelve la columna "nivel" de la tabla de diferencias divididas. Así, para nivel=0,  obtenemos los valores de $y$, si tomamos nivel=1 devuelve las pendientes y así sucesivamente.

#auto def diferencias_divididas_puntos(puntos,nivel): if nivel==0: return [(y,x,x) for (x,y) in puntos] else: dd=diferencias_divididas_puntos(puntos,nivel-1) return [((dd[i+1][0]-dd[i][0])/(dd[i+1][2]-dd[i][1]),dd[i][1],dd[i+1][2]) for i in [0..len(dd)-2]] def diferencias_divididas(puntos,nivel): return [x for x,y,z in diferencias_divididas_puntos(puntos,nivel)] 
       

Para calcular el polinomio de interpolación, recordemos que basta tomar los valores superiores de cada columna, multiplicados por el término correspondiente y sumados. Vamos a calcular el polinomio interpolador de (-1,3), (1,2), (3,-1).

puntos=[(-1,3),(1,2),(3,-1)] p(x)=3+diferencias_divididas(puntos,1)[0]*(x+1)+diferencias_divididas(puntos,2)[0]*(x+1)*(x-1) 
       
points([(-1,3),(1,2),(3,-1)])+plot(p(x),(x,-1,3)) 
       

Practica

Obtén el polinomio de interpolación de los puntos $(1,3)$, $(2,4)$, $(4,1), (5,2)$ usando los dos métodos anteriores.

 
       

Interpolación con funciones de base radial

Interpolación de funciones de una variable

Comenzamos definiendo distintas familias de base radial.

def f_base(r,tipo='gaussiana',epsilon=1.0): if tipo == 'gaussiana': return exp(-(epsilon*r)^2) if tipo == 'distancia': return r if tipo == 'multicuadratica': return sqrt(1+(epsilon*r)^2) if tipo == 'multicuadratica_inversa': return 1/sqrt(1+(epsilon*r)^2) 
       

Para obtener la función de interpolación, tenemos que obtener el valor de los coeficientes. Esto lo hacemos resolviendo el sistema de ecuaciones.

def interpolante(puntos,tipo='gaussiana'): px,py=zip(*puntos) M=Matrix(RDF,[[f_base(abs(x1-x2),tipo=tipo) for x1 in px] for x2 in px]) return sum(map(prod,zip(M.solve_right(py),[f_base(abs(x-cx),tipo=tipo) for cx in px]))) 
       

Vamos a comparar distintas familias.

points(puntos,color='black',size=30)+plot(interpolante(puntos),(x,-2,5),legend_label='gaussiana')+plot(interpolante(puntos,tipo='multicuadratica'),(x,-2,5),color='green',legend_label=u'multicuadrática')+plot(interpolante(puntos,tipo='multicuadratica_inversa'),(x,-2,5),color='red',legend_label=u'multicuadrática inversa') 
       

Interpolación de funciones de dos variables (superficies)

Podemos repetir el proceso para funciones de dos (o más) variables. Para ello, debemos definir la distancia entre los puntos y después resolver el sistema de ecuaciones para obtener los coeficientes.

def d(p1,p2): # Función que devuelve la distancia entre dos puntos return sqrt( (p1[0]-p2[0])^2 + (p1[1]-p2[1])^2 ) var('x,y') def interpolante(puntos,tipo='gaussiana'): px,py,pz=zip(*puntos) M=Matrix(RDF,[[f_base(d(p1,p2),tipo=tipo) for p1 in zip(px,py)] for p2 in zip(px,py)]) return sum(map(prod,zip(M.solve_right(pz),[f_base(d((x,y),p),tipo=tipo) for p in zip(px,py)]))) 
       
puntos=[(-1,-1,1),(-1,1,2),(1,-1,-1),(1,1,2),(0,0,3)] points(puntos,color='black')+plot3d(interpolante(puntos),(x,-2,2),(y,-2,2)) 
       

Practica

Interpola los puntos $(1,3)$, $(2,4)$, $(4,1), (5,2)$ usando los dos métodos anteriores.

 
       

Interpolación mediante splines


Splines lineales

Los splines lineales de una serie de puntos consisten en unirlos mediante una línea recta. Vamos a definir una función auxiliar que nos de la ecuación de la recta que pasa por dos puntos.

#auto def segmento(puntoA,puntoB): return puntoA[1]+(x-puntoA[0])*(puntoB[1]-puntoA[1])/(puntoB[0]-puntoA[0]) 
       

Ahora basta considerar cada intervalo entre dos puntos y calcular el segmento que los une.

puntos=[(-1,3),(1,2),(3,-1)] p1(x)=segmento(puntos[0],puntos[1]) p2(x)=segmento(puntos[1],puntos[2]) points(puntos)+plot(p1(x),(x,-1,1))+plot(p2(x),(x,1,3)) 
       
 
       

Splines cúbicos

Para calcular los splines cúbicos, construiremos una función auxiliar que reciba los puntos y el tipo de spline cúbico y devuelva el sistema de ecuaciones que satisfacen los coeficientes de los polinomios.

#auto def matriz_spline(puntos,tipo='natural'): # Comenzamos creando matrices vacías y definiendo algunos valores MP=Matrix(puntos) n=MP.dimensions()[0]-1 A=matrix(RDF,4*n) b=vector(RDF,4*n) # Rellenamos las condiciones de los puntos extremos A[0,0:4]=matrix(RDF,[[1,MP[0,0],MP[0,0]^2,MP[0,0]^3]]); b[0]=MP[0,1] A[1,(4*n-4):(4*n)]=matrix(RDF,[[1,MP[n,0],MP[n,0]^2,MP[n,0]^3]]); b[1]=MP[n,1] # Rellenamos las condiciones que dependen del tipo de método if tipo=='natural': # La derivada segunda en el origen A[2,0:4]=matrix(RDF,[[0,0,2,6*MP[0,0]]]); b[2]=0 # La derivada segunda en el final A[3,(4*n-4):(4*n)]=matrix(RDF,[[0,0,2,6*MP[n,0]]]); b[3]=0 elif tipo=='periodico': A[2,0:4]=matrix(RDF,[[0,1,2*MP[0,0],6*MP[0,0]]]) A[3,0:4]=matrix(RDF,[[0,0,2,6*MP[0,0]]]) b[2]=0 A[2,(4*n-4):(4*n)]=matrix(RDF,[[0,-1,-2*MP[n,0],-3*MP[n,0]^2]]) A[3,(4*n-4):(4*n)]=matrix(RDF,[[0,0,-2,-6*MP[n,0]]]) b[3]=0 # Rellenamos las condiciones de los puntos intermedios for i in range(1,n): A[4*i,4*(i-1):4*i]=matrix(RDF,[[1,MP[i,0],MP[i,0]^2,MP[i,0]^3]]); b[4*i]=MP[i,1] A[4*i+1,4*i:4*(i+1)]=matrix(RDF,[[1,MP[i,0],MP[i,0]^2,MP[i,0]^3]]); b[4*i+1]=MP[i,1] A[4*i+2,4*(i-1):4*(i+1)]=matrix(RDF,[[0,1,2*MP[i,0],3*MP[i,0]^2, 0,-1,-2*MP[i,0],-3*MP[i,0]^2]]); b[4*i+2]=0 A[4*i+3,4*(i-1):4*(i+1)]=matrix(RDF,[[0,0,2,6*MP[i,0],0,0,-2,-6*MP[i,0]]]); b[4*i+3]=0 return A,b 
       
puntos=zip([0,1,2,3],[4,3,2,4]) A,b=matriz_spline(zip([0,1,2,3],[4,3,2,4])) show(A),show(b) 
       
sol=A.solve_right(b) p1(x)=sol[0]+sol[1]*x+sol[2]*x^2+sol[3]*x^3 p2(x)=sol[4]+sol[5]*x+sol[6]*x^2+sol[7]*x^3 p3(x)=sol[8]+sol[9]*x+sol[10]*x^2+sol[11]*x^3 show(p1),show(p2),show(p3) 
       
points(puntos,size=40)+plot(p1(x),(x,0,1))+plot(p2(x),(x,1,2))+plot(p3(x),(x,2,3)) 
       

Sage dispone de una función para calcular splines, denominada spline

s=spline(puntos) ss=plot(s,(x,0,3),color="blue") cp=points(puntos,size=40,color="blue") show(ss+cp,figsize=4) 
       

Practica

Obtén spline lineal y el spline cúbico natural de los puntos $(1,3)$, $(2,4)$, $(4,1), (5,2)$.

 
       

Interpolación de curvas

Interpolación polinomial de curvas

En el caso de que tengamos una curva definida de forma paramétrica, el procedimiento para interpolar es el siguiente.

Supongamos que tenemos los puntos:

$t=0 \to (1,0)$, $t=1 \to (0,1)$, $t=2 \to (-1,0)$, $t=3\to (0,-1)$, $t=4 \to (0,0)$.

Podemos definir los vectores para cada coordenada:

pt=vector([0,1,2,3,4]);px=vector([1,0,-1,0,0]); py=vector([0,1,0,-1,0]);t=var('t'); 
       

1. Interpolamos los datos del vector (pt,px). En lugar de usar los métodos vistos en el primer apartado, utilizaremos un método de Sage que calcula directamente el polinomio interpolador a partir de los puntos.

R = PolynomialRing(CDF, 't') cx(t)=R.lagrange_polynomial(zip(pt,px)); expand(cx) 
       

2. Interpolamos los datos de (pt,py)

cy(t)=R.lagrange_polynomial(zip(pt,py)); expand(cy) 
       

Dibujamos la curva

parametric_plot((cx(t),cy(t)),(t,0,4),thickness=3)+point(zip(px,py),rgbcolor='red', pointsize=120) 
       


Ajuste por mínimos cuadrados

Regresión Lineal

Supongamos que tenemos el siguiente conjunto de datos:

puntos=[(1,1),(2,2.5),(3,2.75),(4,4.15),(5,4.95),(6,6.1),(7,8),(8,7.5),(9,9.1),(10,9.95)] points(puntos,color='purple',size=10,figsize=4) 
       

Vamos a adaptar un modelo de regresión lineal a los mismos. Para ello usaremos la función find_fit de Sage, que busca los parámetros adecuados para ajustar un modelo a un conjunto de datos, mediante el método de mínimos cuadrados.

var('a,b') model(x)=a+b*x f1=find_fit(puntos,model,solution_dict=True) print 'La recta de regresion es', f1[a]+f1[b]*x points(puntos,color='purple',size=20,figsize=4) + plot(model(a=f1[a],b=f1[b]),(x,0,10),color='red',thickness=2) 
       


Regresión cuadrática

Veamos cómo funciona la regresión cuadrática:

Ajustamos un modelo cuadrático:

var('a,b,c') model2(x)=a+b*x+c*x^2 f2=find_fit(puntos,model2,solution_dict=True) graf2=points(puntos,color='purple',size=20) +plot(model2(a=f2[a],b=f2[b],c=f2[c]),(x,0,10),color='red',thickness=2) show(graf2,figsize=4) print 'Regresión cuadrática', f2[a]+f2[b]*x+f2[c]*x^2 
       


Otras funciones de ajuste

datos=[(91,39.433942),(81,38.341728),(70,35.025399),(60,31.730170),(50,29.552647),(40,27.835548),(30,26.107849),(20,24.203518),(10,23.202388),(0,22.094660)]; points(datos,size=30,figsize=4) 
       

Si ajustamos con un spline cúbico:

s=spline(datos) ss=plot(s,(x,0,91),color="blue") cp=points(datos,size=40,color="blue") show(ss+cp,figsize=4) 
       

Si ajustamos con una función del tipo $a/(b+c e^{-ax+91})$

var('a b c x') assume(a>0, b>0.0001) 
       
model2(x,a,b,c) = a/(b+c*exp(-a*x/91)) 
       
fit=find_fit(datos, model2, parameters = [a, b, c], variables = [x], solution_dict = True) print 'a= ',fit[a]; print 'b= ',fit[b]; print 'c= ',fit[c]; 
       
graf=cp+plot(model2(a=fit[a],b=fit[b],c=fit[c]),(x,0,91),color='red',thickness=2) show(graf, figsize=4) 
       
show(ss+cp+graf,figsize=4) 
       

Practica

Obtén los parámetros a y b  minimicen el error cuadrático de ajustar los puntos $(1,0), (2,4), (4,5), (5,2)$ con una función del tipo f(x)=a+b*x^2.

 
       

Integración Numérica

Vamos a calcular numéricamente la integral $\int_0^1 1+e^{-x}\sin(4*x)\,dx$.

Método del trapecio

El método del trapecio consiste en calcular la integral como el área del trapecio formado por el eje x y los puntos $(a,f(a))$, $(b,f(b)$. Vamos a dibujar la función y el trapecio anterior.

f(x)=1+e^(-x)*sin(4*x) a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(1,f(1))]) b=point2d(v,size=120,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) show(a+b+c+plot(f(0)+(x-0)*(f(1)-f(0)),(x,0,1),fill='axis', color='green',fillcolor='green',thickness=4)) 
       

El valor aproximado de la integral lo calculamos mediante la fórmula:

I1=(f(0.0)+f(1.0))/2.0;I1 
       


Método de Simpson

Por el método de Simpson, la integral se aproxima por

$$I=\frac{h}{3}\left(f(a)+4 f(x_1)+f(b) \right)$$

Es decir:

f(x)=1+exp(-x)*sin(4*x) I2=(f(0.0)+4*f(0.5)+f(1.0))/6;I2 
       

Recordemos que consiste en interpolar la función en los extremos y el punto medio y calcular la integral de dicho polinomio. El polinomio de interpolación es:

R = PolynomialRing(RDF, 'x') px=vector([0,0.5,1]);py=vector([f(0),f(0.5),f(1)]); puntos=zip(px,py); Q(x) = R.lagrange_polynomial(puntos);expand(Q) 
       

Dibujamos ahora la función, el polinomio a integrar y el área que define la integral.

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(0.5,f(0.5)),(1,f(1))]) b=point(v,size=120,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) d=plot(Q(x),(x,0,1),color='green',fill='axis',fillcolor='green',thickness=4) show(a+b+c+d) 
       

Veamos que la integral obtenida por el método de Simpon según la fórmula $A=\frac{h}{3}(f(a)+4f(a+h)+f(b))$ coincide con la integral del polinomio de interpolación:

Q.integrate(x,0,1).n() 
       

Trapecio compuesto

Vamos a calcular la integral mediante el método del trapecio compuesto dividiendo el intervalo en 5 subintervalos.

f(x)=1+e^(-x)*sin(4*x) a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=([(0,f(0)),(0.2,f(0.2)),(0.4,f(0.4)),(0.6,f(0.6)),(0.8,f(0.8)),(1,f(1))]) b=point2d(v,size=60,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) l1=plot(f(0)+(x-0)*(f(0.2)-f(0))*5,(x,0,0.2),fill='axis', color='green',fillcolor='green',thickness=4) l2=plot(f(0.2)+(x-0.2)*(f(0.4)-f(0.2))*5,(x,0.2,0.4),fill='axis', color='green',fillcolor='green',thickness=4) l3=plot(f(0.4)+(x-0.4)*(f(0.6)-f(0.4))*5,(x,0.4,0.6),fill='axis', color='green',fillcolor='green',thickness=4) l4=plot(f(0.6)+(x-0.6)*(f(0.8)-f(0.6))*5,(x,0.6,0.8),fill='axis', color='green',fillcolor='green',thickness=4) l5=plot(f(0.8)+(x-0.8)*(f(1)-f(0.8))*5,(x,0.8,1),fill='axis', color='green',fillcolor='green',thickness=4) a+b+c+l1+l2+l3+l4+l5 
       

Para calcular la integral, aplicamos la fórmula:

I3=(f(0.0)+2*f(0.2)+2*f(0.4)+2*f(0.6)+2*f(0.8)+f(1.0))/10;I3 
       

Integración con Sage

Podemos calcular la integral exacta utilizando el comando "integrate". Si hacemos f.integrate(x,a,b), nos devuelve el valor (exacto) de la integral de la función $f$ entre $a$ y $b$.

f.integrate(x,0,1).n() 
       

En Sage tenemos el comando numerical_integral que nos permite calcular la integral numérica de una función (devuelve la integral numerica y una estimación del error cometido):

numerical_integral(f,0,1) 
       


Practica

Aproxima la integral $\int_{-1}^1 e^{x^2}\,dx$ usando los métodos anteriores.

 
       

Interpolación en superficies

Vamos a suponer que nuestro terreno está descrito por la siguiente función:

var('x,y') f(x,y)=sin(2*pi*x)*exp(cos(2*pi*(x+y)))+cos(2*pi*y) contour_plot(f(x,y),(x,0,1),(y,0,1),fill=False,labels=True,cmap='jet',label_colors='black',label_inline=True) 
       

Supongamos que hemos medido las cotas en una serie de puntos. Para simularlo, tomaremos puntos al azar.

# La siguiente función recibe una lista de valores y los devuelve escalados entre 0 y 1 def normalizar01(lista): minimo=min(lista) maximo=max(lista) return [(k-minimo)/(maximo-minimo) for k in lista] 
       
n=100 puntos=[(random()*1.2-0.1,random()*1.2-0.1) for _ in [1..n]] valores=[f(k[0],k[1]) for k in puntos] colores=normalizar01(valores) 
       

Representamos los puntos y su cota (amarillo para cotas altas y negro para las bajas).

sum([points(puntos[k],color=(colores[k],colores[k],0)) for k in [0..n-1]]) 
       

Realizamos la interpolación. El proceso es un poco complejo: Hay que generar una malla con los puntos donde vamos a evaluar la interpolación. Después pasamos a la función griddata los puntos donde tenemos las medidas, los valores de dichas medidas, la malla donde vamos a interpolar y el método de interpolación. 

import numpy as np from scipy.interpolate import griddata,Rbf ti = np.linspace(0, 1.0, 100) XI, YI = np.meshgrid(ti, ti) ZI = griddata(puntos, valores, (XI, YI), method='cubic') 
       

Finalmente, usaremos la función contour de la librería pyplot para tranformar la interpolación (que es únicamente una matriz con los valores interpolados) en curvas de nivel.

import matplotlib.pyplot as plt plt.figure() CS = plt.contour(XI, YI, ZI) plt.clabel(CS, inline=1) plt.title('Curvas de nivel') plt.savefig("") plt.close() 
       

La interpolación por splines cúbicos no es la única posible. Por ejemplo, podemos interpolar usando el promedio de los valores observados con pesos dados por una función de base radial:

ti = np.linspace(0, 1.0, 100) XI, YI = np.meshgrid(ti, ti) Xp=np.array(zip(*puntos)[0]) Xp=np.array(zip(*puntos)[1]) Zk=np.array([RDF(k) for k in valores]) rbf = Rbf(Xk, Yk, Zk, function='inverse') ZI = rbf(XI, YI) plt.figure() CS = plt.contour(XI, YI, ZI) plt.clabel(CS, inline=1) plt.title('Curvas de nivel') plt.savefig("") plt.close() 
       

Interpolación linear:

ti = np.linspace(0, 1.0, 100) XI, YI = np.meshgrid(ti, ti) ZI = griddata(puntos, valores, (XI, YI),method='linear') plt.figure() CS = plt.contour(XI, YI, ZI) plt.clabel(CS, inline=1) plt.title('Curvas de nivel') plt.savefig("") plt.close() 
       

Podemos ver también la diferencia de tiempo de cálculo:

%time ZI = griddata(puntos, valores, (XI, YI), method='cubic') 
       
%time rbf = Rbf(Xk, Yk, Zk, function='inverse') ZI = rbf(XI, YI) 
       
%time ZI = griddata(puntos, valores, (XI, YI), method='linear') 
       

Si queremos interpolar con Kriging, tenemos una página que ilustra el método y como hacerlo en python y también librerías que lo tienen implementado, como PyKriging y PyKrige.

Practica:

1. Prueba a cambiar el número de puntos medidos (20, 50, 100, 200, 500, 1000, 2000, 5000) y compara los distintos métodos.

2. Prueba con otra función (ej. sin(x^2)+cos(y^2)).