GDIDP_AMAT_16_Practica3

hace 393 días por etlopez16

Interpolación e Integración

1. Interpolación polinomial

1.1 Cálculo del polinomio interpolador mediante un sistema de ecuaciones lineales

Consideremos, por ejemplo, los nodos $(0,4), (1,3), (2,1), (3,4)$. Los definiremos como dos vectores, uno con las coordenadas $x$ y otro con las $y$ y definiremos una lista con los puntos mediante la orden zip (que transforma dos vectores en una lista de puntos):

px=vector([0,1,2,3]);py=vector([4,3,1,4]); puntos=zip(px,py); puntos 
       

Planteamos la matriz dada por las ecuaciones $p(x_i)=y_i$, donde $p$ es el polinomio interpolador y $(x_i,y_i)$ los puntos a interpolar.

A=matrix(RDF,[[1,0,0^2,0^3],[1,1,1^2,1^3],[1,2,2^2,2^3],[1,3,3^2,3^3]]); A 
       

Ahora resolvemos el sistema $Ax=b$, donde $A$ es la matriz anterior y $b$ es el vector de coordenadas $y$:

s=A.solve_right(py); s 
       

El polinomio interpolador es:

P(x)=s[0]+s[1]*x+s[2]*x^2+s[3]*x^3;P 
       

Podemos dibujar los puntos (mediante la orden "point") y el polinomio en una gráfica:

d1=point(puntos,rgbcolor='red', pointsize=80) d2=plot(P(x),(x,-1,3)) show(d1+d2) 
       

Para calcular directamente el polinomio interpolador, Sage utiliza el método de lagrange. Para aplicarlo es necesario primero escribir  R=PolynomialRing(RDF,'x')

y después el polinomio interpolador se calcula como R.lagrange_polynomial(puntos), donde puntos es la lista de puntos a interpolar tal y como la definimos al principio (con la función zip).

R = PolynomialRing(RDF, 'x') px=vector([0,1,2,3]);py=vector([4,3,1,4]); puntos=zip(px,py); Q(x) = R.lagrange_polynomial(puntos);Q 
       

Para comprobar que este vector coincide con el anterior, utilizamos la orden expand, que calcula los paréntesis en el polinomio.

Q.expand() 
       

Practica

1. Calcula el polinomio interpolador de los puntos $(0,1)$, $(1,0)$, $(2,1)$. Comprueba que es correcto.

 
       

1.2 Método de Lagrange

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

Definimos los polinomios de Lagrange:

L0(x)=prod([ x-px[k] for k in [0..3] if k!=0 ])/prod([px[0]-px[k] for k in [0..3] if k!=0]) L1(x)=prod([ x-px[k] for k in [0..3] if k!=1 ])/prod([px[1]-px[k] for k in [0..3] if k!=1]) L2(x)=prod([ x-px[k] for k in [0..3] if k!=2 ])/prod([px[2]-px[k] for k in [0..3] if k!=2]) L3(x)=prod([ x-px[k] for k in [0..3] if k!=3 ])/prod([px[3]-px[k] for k in [0..3] if k!=3]) (plot(L0(x),(x,0,3),color='red')+plot(L1(x),(x,0,3),color='green')+plot(L2(x),(x,0,3),color='purple')+plot(L3(x),(x,0,3),color='blue')).show(figsize=7) 
       

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

p(x)=py[0]*L0(x)+py[1]*L1(x)+py[2]*L2(x)+py[3]*L3(x);p.expand() 
       

1.3 Método de Newton o diferencias divididas

Podemos calcular las diferencias divididas a partir de su definición:

def dd(puntos): # Función que calcula la tabla de diferencias divididas if len(puntos)==1: # Si hay un único punto, devuelve la "y" return [[puntos[0][1].n()]] else: # Si hay más de un punto, aplica la fórmula dinf=dd(puntos[:-1]) dsup=dd(puntos[1:]) numerador=dsup[0][-1]-dinf[0][-1] nuevo_valor=[numerador.n()/(puntos[-1][0]-puntos[0][0])] # El nuevo valor calculado lo ponemos en la primera fila return [dinf[0]+nuevo_valor]+dsup 
       

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. En la tabla anterior, esos valores estan en las posiciones $(i,i)$.

px=vector([0,1,2,3]) py=vector([4,3,1,4]) puntos=zip(px,py) puntos 
       
DD=dd(puntos) DD 
       
P(x)=DD[0][0]+DD[0][1]*(x-0)+DD[0][2]*(x-0)*(x-1)+DD[0][3]*(x-0)*(x-1)*(x-2) P.expand() 
       

Practica

Obtén el polinomio de interpolación de los puntos $(1,3)$, $(2,4)$, $(4,1)$ usando directamente la función lagrange_polynomial, calculando los polinomios de Lagrange y mediante diferencias divididas.

 
       


2. Interpolación a trozos: Splines.

2.1 Splines lineales

 

 

Para una interpolación lineal, basta definir la función a trozos (usando piecewise) tal que en cada intervalo $(x_i,x_{i+1})$, tomamos la recta que une los puntos $(x_i,y_i)$ y $(x_{i+1},y_{i+1})$, cuya ecuación es: \[y=y_i+\frac{y_{i+1}-y_i}{x_{i+1}-x_i}(x-x_i).\]

px=vector([-1,1,2,3]) py=vector([4,3,1,4]) puntos=zip(px,py) splinel=Piecewise([[(px[0],px[1]),py[0]+(x-px[0])*(py[1]-py[0])/(px[1]-px[0])],[(px[1],px[2]),py[1]+(x-px[1])*(py[2]-py[1])/(px[2]-px[1])],[(px[2],px[3]),py[2]+(x-px[2])*(py[3]-py[2])/(px[3]-px[2])]]);splinel 
       
a=points(puntos,rgbcolor='red',pointsize=120) b=splinel.plot(thickness=2) show(a+b,figsize=4) 
       

También se puede calcular sin usar el comando piecewise, guardando la lista de rectas y los intervalos entre los que están definidas.

def r2(p1,p2): # Recta que pasa por p1 y p2 return p1[1]+(x-p1[0])*(p2[1]-p1[1])/(p2[0]-p1[0]) sp_lineal=[(px[k-1],px[k],r2(puntos[k-1],puntos[k])) for k in [1..len(puntos)-1]] sp_lineal 
       
sum([plot(l,(x,a,b)) for a,b,l in sp_lineal]) 
       
 
       

2.2 Splines cúbicos

Veamos cómo calcular el spline cúbico que interpola los nodos $(0,4), (1,3), (2,1), (3,4)$. Sage tiene el comando spline, que devuelve el spline cúbico natural.

px=vector([0,1,2,3]);py=vector([4,3,1,4]); puntos=zip(px,py); puntos S1 = spline(puntos); d1=points(puntos,rgbcolor='red',pointsize=120) d2=plot(S1,(x,0,3)); show(d1+d2) 
       

Desafortunadamente, la orden spline no permite mucho control sobre los resultados (no permite obtener las ecuaciones, por ejemplo).

Para calcular explícitamente la fórmula del spline, podemos formular el sistema de ecuaciones que define el spline cúbico que interpola los valores $(0,4), (1,3), (2,1)$

A=Matrix(QQ,[[1, 0, 0, 0,0,0,0,0],[1 , 1, 1, 1,0,0,0,0],[0 , 0, 0, 0,1,1,1,1],[0 , 0, 0, 0,1,2,4,8],[0 , 1, 2, 3,0,-1,-2,-3],[0 , 0, 2, 6,0,0,-2,-6],[0 , 0, 2, 0,0,0,0,0],[0 , 0, 0, 0,0,0,2,12]]); b=Matrix(QQ,[[4],[3],[3],[1],[0],[0],[0],[0]]);A,b 
       
s1=A.solve_right(b).column(0);s1 
       
p1(x)=s1[0]+s1[1]*x+s1[2]*x^2+s1[3]*x^3;p1 
       
p2(x)=s1[4]+s1[5]*x+s1[6]*x^2+s1[7]*x^3;p2 
       
plot(p1(x),(x,0,1))+plot(p2(x),(x,1,2))+point([(0,4),(1,3),(2,1)],color='red', pointsize=80) 
       

Practica: Calcula el spline lineal que interpola los puntos (1,2), (3,1), (5,7). Calcula también el spline cúbico natural que interpola los puntos. Calcula su valor en $4$ y en $\pi$. ¿Coincide con el lineal?

 
       

3. Interpolación de una función

En caso de que los nodos a interpolar provengan de una función, por ejemplo

\[f(x)=\frac{-4\sqrt{3}}{3}\sin(\frac{\pi x}{3})+\frac{2\sqrt{3}}{3}\sin(\frac{2\pi x}{3})+4,\]

en los puntos $x=0,1,2,3$. Vamos a definir la función y dibujar su gráfica:

a=-4*sqrt(3)/3;b=2*sqrt(3)/3; f(x)=a*sin(pi*x/3)+b*sin(2*pi*x/3)+4 
       
px=vector([0.0,1.0,2.0,3.0]) 
       

El polinomio interpolador debe pasar por los puntos $(x_i,f(x_i))$, que en este caso son:

py=vector([f(0),f(1),f(2),f(3)]); py 
       
(4, 3, 1, 4)
(4, 3, 1, 4)
px=vector([0,1,2,3]);py=vector([f(0),f(1),f(2),f(3)]); puntos=zip(px,py); A=matrix(RDF,[[1,0,0,0],[1,1,1^2,1^3],[1,2,2^2,2^3],[1,3,3^2,3^3]]); s=A.solve_right(py); P(x)=s[0]+s[1]*x+s[2]*x^2+s[3]*x^3;P 
       
x |--> 1.0*x^3 - 3.5*x^2 + 1.5*x + 4.0
x |--> 1.0*x^3 - 3.5*x^2 + 1.5*x + 4.0
p=plot(f,(x,0,3),figsize=4)+plot(P,(x,0,3),color='green')+points(puntos,color='red') show(p) 
       

 

También podemos interpolar una función mediante splines, usando el compando spline:

s = spline(puntos) a=plot(s,(0,3),color='orange') show(a+p) 
       

Gráficamente, es claro que el ajuste proporciona muy bajos errores.