MUI_IIMA_1415_Numerico_3

1047 days ago by etlopez

Interpolació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 
       
[(0, 4), (1, 3), (2, 1), (3, 4)]
[(0, 4), (1, 3), (2, 1), (3, 4)]

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 
       
[ 1.0  0.0  0.0  0.0]
[ 1.0  1.0  1.0  1.0]
[ 1.0  2.0  4.0  8.0]
[ 1.0  3.0  9.0 27.0]
[ 1.0  0.0  0.0  0.0]
[ 1.0  1.0  1.0  1.0]
[ 1.0  2.0  4.0  8.0]
[ 1.0  3.0  9.0 27.0]

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 
       
(4.0, 1.5, -3.5, 1.0)
(4.0, 1.5, -3.5, 1.0)

El polinomio interpolador es:

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

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,0,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) 
       
(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

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

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([4,3,1,4]); puntos=zip(px,py); Q(x) = R.lagrange_polynomial(puntos);Q 
       
x |--> ((x - 3.5)*x + 1.5)*x + 4.0
x |--> ((x - 3.5)*x + 1.5)*x + 4.0

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

Q.expand() 
       
x |--> x^3 - 3.5*x^2 + 1.5*x + 4.0
x |--> x^3 - 3.5*x^2 + 1.5*x + 4.0

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() 
       
x |--> x^3 - 7/2*x^2 + 3/2*x + 4
x |--> x^3 - 7/2*x^2 + 3/2*x + 4

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() 
       
x |--> x^3 - 7/2*x^2 + 3/2*x + 4
x |--> x^3 - 7/2*x^2 + 3/2*x + 4

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 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 
       
Piecewise defined function with 3 parts, [[(-1, 1), -1/2*x + 7/2], [(1,
2), -2*x + 5], [(2, 3), 3*x - 5]]
Piecewise defined function with 3 parts, [[(-1, 1), -1/2*x + 7/2], [(1, 2), -2*x + 5], [(2, 3), 3*x - 5]]
show(splinel) 
       
a=points(puntos,rgbcolor='red',pointsize=120) b=splinel.plot(thickness=4) show(a+b) 
       
 
       

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

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 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]  [4]
[ 1  1  1  1  0  0  0  0]  [3]
[ 0  0  0  0  1  1  1  1]  [3]
[ 0  0  0  0  1  2  4  8]  [1]
[ 0  1  2  3  0 -1 -2 -3]  [0]
[ 0  0  2  6  0  0 -2 -6]  [0]
[ 0  0  2  0  0  0  0  0]  [0]
[ 0  0  0  0  0  0  2 12], [0]
)
(
[ 1  0  0  0  0  0  0  0]  [4]
[ 1  1  1  1  0  0  0  0]  [3]
[ 0  0  0  0  1  1  1  1]  [3]
[ 0  0  0  0  1  2  4  8]  [1]
[ 0  1  2  3  0 -1 -2 -3]  [0]
[ 0  0  2  6  0  0 -2 -6]  [0]
[ 0  0  2  0  0  0  0  0]  [0]
[ 0  0  0  0  0  0  2 12], [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) 
       
 
       

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.

R = PolynomialRing(RDF, 'x') 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) 
       

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


 
       

Ejercicios

$a,b,c, d$ son las cuatro primeras cifras de tu DNI .

1. Calcula el polinomio de interpolación de los puntos $(-1,a),(1,b),(3,c),(5,d)$ 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),(1,b),(3,c),(5,d)$

4. Pinta en una misma gráfica los puntos anteriores, el polinomio de interpolación, el spline lineal y el  spline cúbico natural.