GDIDP_AMAT_1415_Practica4

1052 days ago by etlopez

Interpolación de curvas

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)