MUI_IIMA_1415_Numerico4

1099 days ago by etlopez

Interpolación de curvas y Regresión

1. 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)

R = PolynomialRing(CDF, 't') cx(t)=R.lagrange_polynomial(zip(pt,px)); expand(cx) 
       
t |--> -0.208333333333*t^4 + 1.58333333333*t^3 - 3.29166666667*t^2 +
0.916666666667*t + 1.0
t |--> -0.208333333333*t^4 + 1.58333333333*t^3 - 3.29166666667*t^2 + 0.916666666667*t + 1.0

2. Interpolamos los datos de (pt,py)

cy(t)=R.lagrange_polynomial(zip(pt,py)); expand(cy) 
       
t |--> 0.333333333333*t^3 - 2.0*t^2 + 2.66666666667*t
t |--> 0.333333333333*t^3 - 2.0*t^2 + 2.66666666667*t

Dibujamos la curva

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

2. Interpolación de curvas mediante Splines

Las siguientes funciones permiten calcular splines en Sage y representar curvas a trozos. Para utilizarlas es necesario ejecutar el cuadrado.

La función spline3 recibe una lista de puntos (como las que hemos utilizado en los ejemplos anteriores) y, opcionalmente, el método para calcular el spline y devuelve el spline.

La función curvaTrozos(sx,sy) recibe dos funciones sx y sy (deben ser funciones de la variable x) y dibuja la gráfica dada por (sx(x),sy(x)).

%auto # Función para calcular los splines. Recibe un vector de puntos y devuelve el spline que interpola los puntos # Tiene un parametro metodo que puede tomar los valores lineal, natural (cúbico) o periodico (cúbico) def spline3(puntos,metodo='natural'): if metodo=='lineal': MP=Matrix(puntos) n=MP.dimensions()[0] s=Piecewise([([(MP[k,0],MP[k+1,0]),MP[k,1]+(x-MP[k,0])*(MP[k+1,1]-MP[k,1])/(MP[k+1,0]-MP[k,0])]) for k in range(n-1) ]) return s else: 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] if metodo=='natural': A[2,0:4]=matrix(RDF,[[0,0,2,6*MP[0,0]]]) b[2]=0 A[3,(4*n-4):(4*n)]=matrix(RDF,[[0,0,2,6*MP[n,0]]]) b[3]=0 elif metodo=='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 else: show('No existe ese método.') # Rellenamos las condiciones intermedias 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 # Resolvemos el sistema c=A.solve_right(b); s=Piecewise([([(MP[k,0],MP[k+1,0]),c[4*k]+c[4*k+1]*x+c[4*k+2]*x^2+c[4*k+3]*x^3]) for k in range(n) ]) return s # Función para dibujar una curva dada por dos funciones a trozos def curvaTrozos(sx,sy): Mx=matrix(sx.intervals()) My=matrix(sy.intervals()) n=Mx.dimensions()[0] G=points((sx.functions()[0](x=Mx[0,0]),sy.functions()[0](x=My[0,0])),rgbcolor='red', pointsize=120) for i in range(0,n): G=G+points((sx.functions()[i](x=Mx[i,0]),sy.functions()[i](x=My[i,0])),rgbcolor='red', pointsize=120) G=G+parametric_plot((sx.functions()[i],sy.functions()[i]),(x,sx.intervals()[i][0],sx.intervals()[i][1]),thickness=4) G=G+points((sx.functions()[n-1](x=Mx[n-1,1]),sy.functions()[n-1](x=My[n-1,1])),rgbcolor='red', pointsize=120) return G 
       

Vamos a ilustrar su uso con el ejemplo anterior. Volvemos a definir los vectores de los puntos a interpolar.

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

Interpolamos pt y px

sx=spline3(zip(pt,px));show(sx) 
       

Interpolamos pt y py

sy=spline3(zip(pt,py));show(sy) 
       

Dibujamos la curva interpolada:

show(curvaTrozos(sx,sy)) 
       

3. Curvas de Bezier

3.1 Bezier con un punto de control

Primero vamos a calcular una curva de Bezier con puntos fijos $(0,0)$ y $(1,1)$ y punto de control $(0,1)$.

px=vector(RDF,[0,0,1]);py=vector(RDF,[0,1,1]); puntos=Matrix(zip(px,py));puntos 
       
[0.0 0.0]
[0.0 1.0]
[1.0 1.0]
[0.0 0.0]
[0.0 1.0]
[1.0 1.0]

Sage admite que escribamos diréctamente la fórmula y nos devuelve la curva como una función de dos componentes:

B=(1-x)^2*puntos[0]+2*(1-x)*x*puntos[1]+x^2*puntos[2]; B 
       
(x^2, -2.0*(x - 1)*x + x^2)
(x^2, -2.0*(x - 1)*x + x^2)

Dibujamos ahora la curva $B$ (en azul), los nodos y puntos de control en rojo y las líneas que unen los nodos con el punto de control (negro). Podemos observar que la curva está dentro del triángulo definido por los tres puntos y que es tangente a las líneas que unen los nodos con el punto de control.

point(zip(px,py),color='red', pointsize=120)+parametric_plot(B,(x,0,1),thickness=4)+line([puntos[0],puntos[1]],color='black',thickness=3)+line([puntos[1],puntos[2]],color='black',thickness=3) 
       

3.2 Bezier con dos puntos de control

Vamos ahora a considerar la curva de Bezier con dos nodos y dos puntos de control.

px=vector(RDF,[0,1,2,3]);py=vector(RDF,[4,3,1,4]); puntos=Matrix(zip(px,py));puntos 
       
[0.0 4.0]
[1.0 3.0]
[2.0 1.0]
[3.0 4.0]
[0.0 4.0]
[1.0 3.0]
[2.0 1.0]
[3.0 4.0]

Usamos la fórmula para cuatro puntos.

B=(1-x)^3*puntos[0]+3*(1-x)^2*x*puntos[1]+3*(1-x)*x^2*puntos[2]+x^3*puntos[3]; B 
       
(3.0*(x - 1)^2*x - 6.0*(x - 1)*x^2 + 3.0*x^3, -4.0*(x - 1)^3 + 9.0*(x -
1)^2*x - 3.0*(x - 1)*x^2 + 4.0*x^3)
(3.0*(x - 1)^2*x - 6.0*(x - 1)*x^2 + 3.0*x^3, -4.0*(x - 1)^3 + 9.0*(x - 1)^2*x - 3.0*(x - 1)*x^2 + 4.0*x^3)

Y dibujamos de nuevo la curva $B$ (en azul), los nodos y puntos de control en rojo y las líneas que unen los nodos con los puntos de control (negro). Podemos observar que la curva está dentro del cuadrilátero definido por los cuatro puntos y que es tangente a las líneas que unen los nodos con los puntos de control.

point(zip(px,py),color='red', pointsize=120)+parametric_plot(B,(x,0,1),thickness=4)+line([puntos[0],puntos[1]],color='black',thickness=3)+line([puntos[2],puntos[3]],color='black',thickness=3) 
       

Un segundo ejemplo

Vamos a ilustrar la flexibilidad de las curvas de Bezier con un segundo ejemplo.

px=vector(RDF,[0,2,-1,1]);py=vector(RDF,[0,2,2,0]); puntos=Matrix(zip(px,py)); B=(1-x)^3*puntos[0]+3*(1-x)^2*x*puntos[1]+3*(1-x)*x^2*puntos[2]+x^3*puntos[3]; B point(zip(px,py),color='red',pointsize=120)+parametric_plot(B,(x,0,1),thickness=4)+line([puntos[0],puntos[1]],color='black',thickness=3)+line([puntos[2],puntos[3]],color='black',thickness=3) 
       
 
       

Ejercicio 1. Obten la interpolación polinomial de la curva que pasa por los puntos $(0,1)$, $(1,1)$, $(1,0)$, $(0,0)$ en los tiempos $0,1,2,3$ y dibújala.

 
       

Ejercicio 2. Calcula la curva de Bezier dada por los nodos $(0,0)$ y $(0,0)$ (el mismo punto repetido) y puntos de control $(1,1)$ y $(-1,1)$. ¿Se puede conseguir que la curva obtenida sea cerrada y diferenciable (es decir, que no tenga "picos")? Razona la respuesta.

 
       

Problema 1.

Interpola la curva paramétrica definida a continuación mediante splines. Para ello elige al menos 40 puntos para interpolar y aplica el método anterior. Elige el tipo de splines que consideres más adecuado para esta interpolación. 

#Curva a interpolar cx(t)=sin(t)*(e^(cos(t))-2*cos(4*t)-sin(t/12)^5) cy(t)=cos(t)*(e^(cos(t))-2*cos(4*t)-sin(t/12)^5) parametric_plot((cx(t),cy(t)),(t,-pi,pi)) 
       

4. Regresión

4.1 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) 
       
La recta de regresion es 0.9890909090908853*x + 0.15999999999816839
La recta de regresion es 0.9890909090908853*x + 0.15999999999816839

Análisis del error:  

residuos=vector(RDF,10); sumares=0; for i in [0..9]: residuos[i]=f1[a]+f1[b]*puntos[i][0]-puntos[i][1]; sumares=sumares+residuos[i]^2 print 'Residuos: ', residuos; print 'Suma de errores al cuadrado: ',sumares 
       
Residuos:  (0.149090909089, -0.36181818182, 0.377272727271,
-0.0336363636383, 0.155454545453, -0.00545454545652, -0.916363636366,
0.572727272725, -0.0381818181839, 0.100909090907)
Suma de errores al cuadrado:  1.50018181818
Residuos:  (0.149090909089, -0.36181818182, 0.377272727271, -0.0336363636383, 0.155454545453, -0.00545454545652, -0.916363636366, 0.572727272725, -0.0381818181839, 0.100909090907)
Suma de errores al cuadrado:  1.50018181818

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 
       

Regresión cuadrática -0.0075757575779553665*x^2 + 1.0724242424244026*x -
0.006666666668869725

Regresión cuadrática -0.0075757575779553665*x^2 + 1.0724242424244026*x - 0.006666666668869725

Análisis del error:

residuos=vector(RDF,10); sumares=0; for i in [0..9]: residuos[i]=f2[a]+f2[b]*puntos[i][0]+f2[c]*(puntos[i][0])^2-puntos[i][1]; sumares=sumares+residuos[i]^2 print 'Residuos o errores: ', residuos; print 'Suma de errores al cuadrado: ',sumares 
       
Residuos o errores:  (0.0581818181776, -0.392121212132, 0.392424242403,
0.0118181817815, 0.216060606004, 0.0551515150712, -0.870909091018,
0.587878787737, -0.0684848486636, 0.00999999977962)
Suma de errores al cuadrado:  1.46987878788
Residuos o errores:  (0.0581818181776, -0.392121212132, 0.392424242403, 0.0118181817815, 0.216060606004, 0.0551515150712, -0.870909091018, 0.587878787737, -0.0684848486636, 0.00999999977962)
Suma de errores al cuadrado:  1.46987878788

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]; 
       
a=  0.397279745165
b=  -0.00742887280603
c=  0.0257466813518
a=  0.397279745165
b=  -0.00742887280603
c=  0.0257466813518
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) 
       

Problema: Ajustar mediante regresión el conjunto de datos correspondiente al número de horas de sol que tiene cada día del año en Badajoz, medidos en la estación meteorológica de Talavera, a una función del tipo

$$\mathbf{a\cos(2bx\pi/365+c)+d}$$

1. Pintar los datos que se proporcionan a continuación (vector de puntos: cord)

2. Definir el modelo para la función de ajuste $a*cos(b*x*2*pi/365+c)+d$

3. Ajustar los datos al modelo y reflejar el ajuste en una gráfica

4. Utilizar el modelo para predecir cuántas horas de sol habrá el próximo 1 de diciembre.

Horas de sol de cada día del año

var('a b c d x') import scipy import numpy time=[9.5, 9.5, 9.51667, 9.53333, 9.55, 9.56667, 9.58333, 9.6, 9.61667, 9.63333, 9.65, 9.66667, 9.68333, 9.71667, 9.73333, 9.75, 9.8, 9.81667, 9.83333, 9.86667, 9.88333, 9.91667, 9.95, 9.96667, 10, 10.05, 10.0667, 10.1, 10.1333, 10.1667, 10.1833, 10.2333, 10.2667, 10.3, 10.3333, 10.3667, 10.4, 10.45, 10.4833, 10.5167, 10.55, 10.5833, 10.6167, 10.6667, 10.7167, 10.75, 10.7833, 10.8167, 10.8667, 10.9, 10.9333, 10.9833, 11.0167, 11.05, 11.1167, 11.15, 11.1833, 11.2333, 11.2667, 11.3167, 11.35, 11.4, 11.4333, 11.4667, 11.5167, 11.55, 11.6, 11.6333, 11.6833, 11.7333, 11.7667, 11.8167, 11.85, 11.9, 11.9333, 11.9833, 12.0167, 12.0667, 12.1167, 12.15, 12.2333, 12.2833, 12.3167, 12.3667, 12.4, 12.45, 12.5, 12.5333, 12.5833, 12.6167, 12.6667, 12.7, 12.75, 12.8, 12.8333, 12.8833, 12.9167, 12.95, 12.9833, 13.0333, 13.0667, 13.1167, 13.15, 13.2, 13.2333, 13.2667, 13.3167, 13.35, 13.4, 13.4333, 13.4667, 13.5167, 13.55, 13.5833, 13.6333, 13.6667, 13.7, 13.75, 13.7833, 13.8167, 13.85, 13.8833, 13.9333, 13.9667, 14, 14.0167, 14.05, 14.0833, 14.1167, 14.15, 14.1833, 14.2167, 14.25, 14.2833, 14.3167, 14.35, 14.3833, 14.4, 14.4167, 14.45, 14.4833, 14.5, 14.5333, 14.5667, 14.5667, 14.6, 14.6167, 14.65, 14.65, 14.6833, 14.7, 14.7333, 14.7333, 14.75, 14.7833, 14.7833, 14.8, 14.8, 14.8333, 14.8333, 14.85, 14.85, 14.8667, 14.8667, 14.8833, 14.8833, 14.8833, 14.9, 14.9, 14.9, 14.8833, 14.9, 14.9, 14.9, 14.8833, 14.8833, 14.8833, 14.8667, 14.8667, 14.8667, 14.85, 14.85, 14.85, 14.8333, 14.8167, 14.7833, 14.7833, 14.7667, 14.75, 14.7333, 14.7167, 14.7, 14.6833, 14.65, 14.65, 14.6167, 14.6, 14.5667, 14.55, 14.5333, 14.5, 14.4667, 14.45, 14.4333, 14.4, 14.3667, 14.3333, 14.3167, 14.2833, 14.25, 14.2333, 14.2, 14.1667, 14.1333, 14.1, 14.05, 14.0167, 13.9833, 13.9667, 13.9333, 13.9, 13.8667, 13.8167, 13.7833, 13.75, 13.7, 13.6667, 13.6333, 13.6, 13.5667, 13.5333, 13.4833, 13.45, 13.4167, 13.3667, 13.3333, 13.2833, 13.25, 13.2, 13.1833, 13.1333, 13.1, 13.05, 13.0167, 12.9667, 12.9333, 12.8833, 12.85, 12.8, 12.7667, 12.7333, 12.6833, 12.65, 12.6, 12.5667, 12.5167, 12.4667, 12.4333, 12.4, 12.3667, 12.3167, 12.2667, 12.2333, 12.1833, 12.1333, 12.1, 12.05, 12.0167, 11.9833, 11.9333, 11.9, 11.85, 11.8167, 11.7667, 11.7333, 11.6833, 11.6333, 11.6, 11.55, 11.5167, 11.4667, 11.4667, 11.4333, 11.3833, 11.35, 11.3, 11.2667, 11.2167, 11.1833, 11.15, 11.1, 11.0667, 11.0167, 10.9833, 10.95, 10.9, 10.8667, 10.8333, 10.7833, 10.75, 10.7167, 10.6833, 10.6333, 10.6, 10.5667, 10.5333, 10.5, 10.45, 10.4167, 10.3833, 10.35, 10.3167, 10.2833, 10.25, 10.2167, 10.1833, 10.15, 10.1, 10.0833, 10.05, 10.0167, 9.98333, 9.96667, 9.93333, 9.91667, 9.88333, 9.85, 9.81667, 9.8, 9.76667, 9.75, 9.71667, 9.7, 9.68333, 9.65, 9.63333, 9.61667, 9.6, 9.58333, 9.58333, 9.56667, 9.55, 9.53333, 9.51667, 9.5, 9.48333, 9.48333, 9.46667, 9.46667, 9.46667, 9.45, 9.45, 9.45, 9.43333, 9.45, 9.43333, 9.45, 9.43333, 9.45, 9.45, 9.46667, 9.46667, 9.46667, 9.46667, 9.48333, 9.5] cord=[] for i in range(365): cord.append((i+1,time[i]))