GDIDP_AMAT_1415_Practica3

1115 days ago by etlopez

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) 
       

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.

 
       

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 
       
Traceback (click to the left of this block for traceback)
...
NameError: name 'puntos' is not defined
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_5.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("UiA9IFBvbHlub21pYWxSaW5nKFJERiwgJ3gnKSAjIFAgZXMgdW4gcG9saW5vbWlvCkREPVIuZGl2aWRlZF9kaWZmZXJlbmNlKHB1bnRvcywgZnVsbF90YWJsZT1UcnVlKTtzaG93KEREKSAjIFRhYmxhIGRlIGRpZmVyZW5jaWFzIGRpdmlkaWRhcw=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpVvBNDZ/___code___.py", line 3, in <module>
    exec compile(u'DD=R.divided_difference(puntos, full_table=True);show(DD) # Tabla de diferencias divididas
  File "", line 1, in <module>
    
NameError: name 'puntos' is not defined

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 directamente la función lagrange_polynomial, calculando los polinomios de Lagrange y mediante diferencias divididas.

 
       

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=2) 
       

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) 
       
x |--> -0.961538461538*x^2 + 1.0
x |--> -0.961538461538*x^2 + 1.0
b=plot(p2(x),(x,-1,1),color='red',thickness=2) 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) 
       
x |--> 2.40384615385*x^4 + (4.4408920985e-16)*x^3 - 3.36538461538*x^2
+ (3.33066907388e-16)*x + 1.0
x |--> 2.40384615385*x^4 + (4.4408920985e-16)*x^3 - 3.36538461538*x^2 + (3.33066907388e-16)*x + 1.0
c=plot(p4(x),(x,-1,1),color='green',thickness=2) p=point2d(v,size=40,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) 
       
x |--> -30.0480769231*x^6 - (3.5527136788e-15)*x^5 +
43.2692307692*x^4 - (7.1054273576e-15)*x^3 - 14.1826923077*x^2 -
(5.55111512312578e-16)*x + 1.0
x |--> -30.0480769231*x^6 - (3.5527136788e-15)*x^5 + 43.2692307692*x^4 - (7.1054273576e-15)*x^3 - 14.1826923077*x^2 - (5.55111512312578e-16)*x + 1.0
d=plot(p6(x),(x,-1,1),color='orange',thickness=2) 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) 
       
x |--> 260178.630972*x^20 - 1012094.87448*x^18 -
(1.51339918375e-09)*x^17 + 1639177.41085*x^16 - (1.42026692629e-08)*x^15
- 1442944.96344*x^14 - (1.51339918375e-08)*x^13 + 757287.092775*x^12 -
(2.85217538476e-09)*x^11 - 245249.283925*x^10 - (3.63797880709e-11)*x^9
+ 49317.5429129*x^8 - (2.09183781408e-11)*x^7 - 6119.21994942*x^6 -
(1.64845914696e-12)*x^5 + 470.846226759*x^4 + (2.84217094304e-14)*x^3 -
24.1434796248*x^2 - (2.22044604925031e-16)*x + 1.0
x |--> 260178.630972*x^20 - 1012094.87448*x^18 - (1.51339918375e-09)*x^17 + 1639177.41085*x^16 - (1.42026692629e-08)*x^15 - 1442944.96344*x^14 - (1.51339918375e-08)*x^13 + 757287.092775*x^12 - (2.85217538476e-09)*x^11 - 245249.283925*x^10 - (3.63797880709e-11)*x^9 + 49317.5429129*x^8 - (2.09183781408e-11)*x^7 - 6119.21994942*x^6 - (1.64845914696e-12)*x^5 + 470.846226759*x^4 + (2.84217094304e-14)*x^3 - 24.1434796248*x^2 - (2.22044604925031e-16)*x + 1.0
e=plot(p20(x),(x,-1,1),color='violet',thickness=2) 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)$.

 
       

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

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 
       
([ 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], [4]
[3]
[3]
[1]
[0]
[0]
[0]
[0])
([ 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], [4]
[3]
[3]
[1]
[0]
[0]
[0]
[0])
s1=A.solve_right(b).column(0);s1 
       
(4, -3/4, 0, -1/4, 7/2, 3/4, -3/2, 1/4)
(4, -3/4, 0, -1/4, 7/2, 3/4, -3/2, 1/4)
p1(x)=s1[0]+s1[1]*x+s1[2]*x^2+s1[3]*x^3;p1 
       
x |--> -1/4*x^3 - 3/4*x + 4
x |--> -1/4*x^3 - 3/4*x + 4
p2(x)=s1[4]+s1[5]*x+s1[6]*x^2+s1[7]*x^3;p2 
       
x |--> 1/4*x^3 - 3/2*x^2 + 3/4*x + 7/2
x |--> 1/4*x^3 - 3/2*x^2 + 3/4*x + 7/2
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)