GIIT_AMA_1415_Tema_1_Practica_3

1080 days ago by trinidad1415

Interpolación polinomial

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 points) y el polinomio en una gráfica:

d1=points(puntos,rgbcolor='red', pointsize=80) d2=plot(P(x),(x,0,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() 
       

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(RDF,[0,1,2,3]);py=vector(RDF,[f(0),f(1),f(2),f(3)]); puntos=zip(px,py); P(x)=R.lagrange_polynomial(puntos);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 $1$ al $5$. La llamaremos $lexponentes$
  • Escribe [1]+[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.
 
       

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

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.

px=vector([-1,1,2,3]);py=vector([4,3,1,4]); puntos=zip(px,py); 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

 
       

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)