GT_MM_1415_Tema_3_Interpolacion_Ajuste_Integracion

1041 days ago by trinidad1415

Interpolación polinomial

Método de Lagrange

Vamos a calcular el polinomio interpoladorusando 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)*(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 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.