GIIT_AMA_1314_Numerico_Interpolacion_Integracion

1502 days ago by trinidad

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([-1,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,-1,(-1)^2,(-1)^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) 
       

1.2 Interpolación de una función

Supongamos que se trata de interpolar 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 
       

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

f(0),f(1),f(2),f(3) 
       
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 
       

Es decir, debemos obtener el polinomio que interpola los puntos $(0,4),(1,3),(2,1),(3,4)$. Son los mismos puntos que en el apartado anterior y por tanto ya sabemos que el polinomio interpolador es \[p(x)=1.0 \, x^{3} - 3.5 \, x^{2} + 1.5 \, x + 4.0\]

Dibujamos la función y el polinomio interpolador:

d1=point(puntos,rgbcolor='red', pointsize=80) d2=plot(P(x),(x,0,3)) d3=plot(f(x),x,0,3,color='green') (d1+d2+d3).show(figsize=7) 
       

Podemos dibujar también el error cometido. Vemos que para los puntos interpolados (puntos entre datos conocidos) el error es pequeño, pero para los puntos extrapolados (puntos fuera del intervalo donde hay datos) el error crece rápidamente.

plot(f(x)-P(x),(x,-1,4)) 
       

Practica

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

2. Vamos a calcular el polinomio interpolador de la función $\sin$ en el intervalo $[0,\pi/2]$. Para ello tomaremos 6 puntos equiespaciados.

  • Crea un vector con los 6 puntos equiespaciados comenzando en $0$ y terminando en $\pi/2$. Para ello utiliza range con el punto inicial, el final y paso $\pi/10$. Vamos a llamarla $lx$
  •  Utilizando la función map aplicada a la lista anterior, calcula $\sin x_i$ para cada $x_i$ de la lista anterior.
  • Crea ahora una lista con los números del $0$ al $5$. La llamaremos $lexponentes$
  • Escribe [lx[0]^k for k in lexponentes]. Como verás has obtenido la primera fila de la matrix de interpolación.
  • Para crear la matriz entera, repite el paso anterior, o crea una función que lo haga, o escribe la expresión dentro de otros corchetes para que se repita por cada elemento de $lx$ ([ expresión_anterior_modificada for xi in lx]).
  • Como ya tenemos todos los elementos, podemos resolver el sistema. Sea $cl$ la solución
  • Ahora hacemos P(x)=sum([cl[k]*x^k for k in lexponentes]) y tenemos el polinomio.
  • Representamos la función, el polinomio y los puntos con las órdenes anteriores.
 
       

1.3. Método de Lagrange

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([1,2,3,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() 
       

Este proceso también se podría hacer "artesanalmente", definiendo 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.4 Diferencias divididas o método de interpolación de Newton

El comando divided_difference calcula la tabla de diferencias divididas. Devuelve una lista de vectores tal que si tomamos las primeras posiciones de cada vector obtenemos los valores de $y$ (primera columna), si tomamos las segundas posiciones, las diferencias divididas de primer nivel (segunda columna en nuestra tabla de diferencias dividias) y así sucesivamente.

R = PolynomialRing(RDF, 'x') # P es un polinomio DD=R.divided_difference(puntos, full_table=True);show(DD) # Tabla de diferencias divididas 
       

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

P(x)=DD[0][0]+DD[1][1]*(x-0)+DD[2][2]*(x-0)*(x-1)+DD[3][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 diréctamente la función lagrange_polynomial, calculando los polinomios de Lagrange y mediante diferencias divididas.

Para los polinomios de Lagrange, se recomienda generar una función Lagrange(x,i) que devuelva el polinomio $L_i$. Se puede hacer escribiendo:

def Lagrange(x,i):

  return escribir_la_definicion_de_polinomio_de_Lagrange

 
       

1.5 Fenómeno de Runge

En ocasiones, subir el número de puntos y, por tanto, el grado del polinomio interpolador de una función no disminuye el error cometido sino que lo aumenta. Veámoslo con un ejemplo. Vamos a considerar la función \[f(x)=\frac{1}{1+25x^2}\] en el intervalo $(-1,1)$.

f(x)=1/(1+25*x^2) 
       
a=plot(f(x),(x,-1,1),color='blue',thickness=4) 
       

Tomemos tres puntos e interpolemos la función.

v=([(-1,f(-1)),(0,f(0)),(1,f(1))]) p2(x)=R.lagrange_polynomial(v); expand(p2) 
       
b=plot(p2(x),(x,-1,1),color='red',thickness=4) p=point2d(v,size=60,color='black') show(p+a+b+p) 
       

Tomamos ahora 5 puntos.

v=([(-1,f(-1)),(-0.6,f(-0.6)),(0,f(0)),(0.6,f(0.6)),(1,f(1))]) p4(x)=R.lagrange_polynomial(v); expand(p4) 
       
c=plot(p4(x),(x,-1,1),color='green',thickness=4) p=point2d(v,size=60,color='black') show(a+b+c+p) 
       

Podemos observar que el error cometido ha aumentado en algunos puntos. Vamos a añadir más puntos para interpolar.

v=([(-1,f(-1)),(-0.6,f(-0.6)),(-0.2,f(-0.2)),(0,f(0)),(0.2,f(0.2)),(0.6,f(0.6)),(1,f(1))]) p6(x)=R.lagrange_polynomial(v); expand(p6) 
       
d=plot(p6(x),(x,-1,1),color='orange',thickness=4) p=point2d(v,size=60,color='black') show(a+b+c+d+p) 
       

Finalmente, tomemos 21 puntos para hacer la interpolación.

v=([(-1,f(-1)),(-0.9,f(-0.9)),(-0.8,f(-0.8)),(-0.7,f(0.7)),(-0.6,f(-0.6)),(-0.5,f(-0.5)),(-0.4,f(-0.4)),(-0.3,f(0.3)),(-0.2,f(-0.2)),(-0.1,f(0.1)),(0,f(0)),(0.1,f(0.1)),(0.2,f(0.2)),(0.3,f(0.3)),(0.4,f(0.4)),(0.5,f(0.5)),(0.6,f(0.6)),(0.7,f(0.7)),(0.8,f(0.8)),(0.9,f(0.9)),(1,f(1))]) p20(x)=R.lagrange_polynomial(v); expand(p20) 
       
e=plot(p20(x),(x,-1,1),color='violet',thickness=4) p=point2d(v,size=30,color='black') show(a+b+c+d+e+p,ymin=-5,ymax=5) 
       


 
       

2. Interpolación a trozos: Splines.

2.1. 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),...]).


 
       

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([-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:

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) 
       
 
       

4. Integración Numérica

Vamos a calcular numéricamente la integral $\int_0^1 1+e^{-x}\sin(4x)\,dx$.

4.1 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)) 
       

4.1 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() 
       

4.3 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 
       

4.4 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) 
       

Ejercicios

$a,b,c$ son las tres primeras cifras de tu DNI (del DNI medio si sois dos).

1. Calcula el polinomio de interpolación de los puntos $(-1,a),(2,b),(4,c)$ con los tres métodos vistos anteriormente.

2. Calcula el polinomio interpolador de la función $e^{ax^2-bx+c}$, en $11$ puntos equiespaciados del intervalo $[0,1]$.

3. Calcula el spline lineal que interpola $(-1,a),(2,b),(4,c).

4. Calcula la integral

$$\int_0^\pi \frac{a \sin x+b\cos x}{10+c\sin x} dx$$

con los métodos del trapecio, Simpson y trapecio con cinco intervalos. Compara el error cometido con la integral exacta.