GIIT_AMA_1415_Tema_1_Practica_4

1140 days ago by trinidad1415

Interpolación a trozos: Splines.

Splines lineales

Para definir los splines lineales, necesitamos poder definir las funciones a trozos. En Sage, lo hacemos mediante el comando piecewise

 

f = Piecewise([[(-1,1),x^2],[(1,3),1-x]]); plot(f) 
       

Para una interpolación lineal, basta definir la función a trozos 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); puntos 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=4) show(a+b) 
       

Practica

Crea el spline lineal que interpola los puntos $(0,1)$, $(1,0)$ y $(2,1)$.

Para ello, lo más sencillo es crear una función que reciba dos puntos: $(x1,y1)$, $(x2,y2)$ y devuelva el intervalo de $x$ y la recta que pasa por los dos puntos:

 [(x1,x2),y1+(x-x1)*(y2-y1)/(x2-x1)]

La función la puedes definir como

def r2(x1,y1,x2,y2):

   return ...

Una vez tengas la función, puedes crear el spline con Piecewise([r2(x0,y0,x1,y1),r2(x1,y1,x2,y2),...]).


 
       

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([-1,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,-1,3)); show(d1+d2) 
       

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

Veamos cómo calcular el spline cúbico que interpola los nodos $(0,4), (1,3), (2,1)$, para lo cual debemos importar la función spline

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) 
       

 

También podemos interpolar una función, como por ejemplo, la función seno, mediante splines, usando el compando spline:

values = [ (x,sin(x)) for x in range(10)] interpolation = spline(values) a=plot(interpolation,(0,9)) + points(values,color='red', pointsize=80) show(a) 
       

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

g(x)=sin(x) b=plot(g(x),(x,0,9),color='green') show(a+b) 
       

Practica

1. Calcula el spline lineal que interpola los puntos (1,2), (3,1), (5,7). ¿Cuanto vale en 4? ¿y en $\pi$?

2. Calcula el spline cúbico natural que interpola los puntos del ejercicio anterior. Calcula su valor en $4$ y en $\pi$. ¿Coincide con el lineal?

3. (avanzado) Vamos a crear una función que nos calcule la matriz de los splines cúbicos naturales.

  • Creamos la cabecera de la función (def matrz_splines(puntos):)
  • Tomamos n como el número de puntos (n=len(puntos))
  • Creamos una matriz de ceros M (M=Matrix(RDF,4*(len(puntos)-1)))
  • Vamos a hacer cuatro bucles, cada uno guardará una de las condiciones que cumple un spline.
  • El primer bucle guardará en las filas 0 hasta n-2 de M la condición de que interpole el extremo izquierdo. Es conveniente guardar el valor de la x del punto que estamos considerando (vx=puntos[i][0]) y entonces basta seleccionar las cuatro posiciones de M y guardar los cuatro valores (M[i,i*4:i*4+4]=[[1,vx,vx^2,vx^3]]). Es conveniente incluir return M en este punto y comenzar a probar la función.
  • El segundo bucle guardará en las filas n-1 hasta 2*n-3 la condición de que interpole el extremo izquierdo.
  • El tercer bucle guardará en las filas 2*n-2 hasta 3*n-5 de M la condición de que el spline sea diferenciable. Por ejemplo, se puede escribir así:

    for i in [0..n-3]:
        vx=puntos[i+1][0]
        M[i+2*(n-1),i*4:i*4+4]=[[0,1,2*vx,3*vx^2]]
        M[i+2*(n-1),i*4+4:i*4+8]=[[0,-1,-2*vx,-3*vx^2]]

  • El cuarto bucle guardará en las filas 3*n-4 hasta 4*n-7 de M la condición de que el spline
  • Finalmente en las dos últimas filas guardaremos la condición de que sea natural.
  • Devolvemos M
 
       

Integración Numérica

Vamos a calcular numéricamente la integral $\int_0^1 1+e^{-x}\sin(4x)\,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))$. Aplicamos la fórmula:

f(x)=1+exp(-x)*sin(4*x) I1=(f(0.0)+f(1.0))/2.0;I1 
       

Vamos a dibujar la función y el trapecio anterior.

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)) #show(a+b+line([(1,0),(1,f(1))],color='black',thickness=4)+line([(0,f(0)),(1,f(1))],color='red',thickness=4)) 
       

 

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+exp(-x)*sin(4*x) 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 
       

Y representamos la función y el área calculada

a=plot(f(x),(x,0,1),color='blue',thickness=4,ymin=0,ymax=2) v=[(k,f(k)) for k in [0,0.2..1]] b=point2d(v,size=60,color='black') c=line([(1,0),(1,f(1))],color='black',thickness=2) l=sum([plot(f(k)+(x-k)*(f(k+0.2)-f(k))*5,(x,k,k+0.2),fill='axis', color='green',fillcolor='green',thickness=4) for k in [0,0.2.. 0.8]]) a+b+c+l 
       

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

1. Calcula la integral:

$$\int_0^1 e^{x^2}\, dx $$

con los métodos anteriores. Compara los resultados obtenidos.

2. (avanzado) Escribe una función que reciba una función a integrar, un intervalo y $n$ y devuelva la aproximación a la integral de la función en el intervalo por el método del trapecio compuesto con $n$ intervalos. Para generar los puntos en lo que tienes que evaluar, puedes utilizar [a,(a+(b-a)/n)..b]. Un poco más avanzado, que la función tenga exáctamente dos líneas, la primera (def ...) y la segunda (return ...).