Práctica 8. Interpolación
En esta práctica vamos a ver cómo definir y utilizar polinomios en Matlab y cómo hallar con MatLab el polinomio de interpolación de Lagrange, es decir, el polinomio de grado menor o igual que n que interpola n+1 datos de Lagrange. Este polinomio lo hallaremos resolviendo un sistema de ecuaciones lineales y, también, usando la forma de Newton de dicho polinomio. Para esta última necesitaremos aprender a construir con MatLab la tabla de diferencias divididas.
Contents
Polinomios en Matlab
En Matlab se pueden expresar polinomios mediante vectores. Concretamente el polinomio p1(x) = aN*x^N + ... + a2*x^2 + a1*x +a0 se expresa mediante un vector con sus coeficientes, p1 = [aN, ..., a2, a1, a0]. Por ejemplo el polinomio p(x) = x^3 + 3*x^2 - 7*x + 5 se expresa como
p=[1 3 -7 5]
p =
1 3 -7 5
Para evaluar un polinomio en x0 empleamos el comando polyval( polinomio, x0), donde x0 puede ser un valor numérico o un vector. Por ejemplo para evaluar el polinomio p en 1 y obtener p(1)=1^3 + 3*1^2 - 7*1 + 5 = 2 usamos la instrucción
polyval(p,1)
ans =
2
Si quisieramos representar gráficamente la función p(x) podríamos hacerlo definiendo un vector de abscisas xD y evaluando p en ese vector para obtener un vector de ordenadas yD, del siguiente modo:
xD=[-3:0.1:3]; yD=polyval(p,xD); plot(xD,yD)
Para hallar las raíces de un polinomio p podemos emplear el comando roots. La siguiente instrucción obtiene las raíces de p, vemos que una de ellas es real y las otras dos son complejas.
roots(p)
ans = -4.7111 + 0.0000i 0.8556 + 0.5739i 0.8556 - 0.5739i
Para sumar dos polinomios debemos tener en cuenta sus grados. Debemos rellenar con ceros el polinomio de menor grado, y así sumar dos vectores de la misma dimensión. Veamos como efectuar la siguiente suma : (x^5 + 2*x^3 + x) + (x^2 + 2*x +1) = x^5 + 2*x^3 + x^2 + 3*x + 1 .
p1=[1, 0, 2, 0, 1, 0]; p2=[1, 2, 1]; p1+[ 0, 0, 0, p2]
ans =
1 0 2 1 3 1
Para multiplicar dos polinomios emplearemos el comando conv. Así para efectuar la multiplicación (x^2 + 2*x + 3)*(-x^3 + 1) = -x^5 - 2*x^4 - 3*x^3 + x^2 + 2*x + 3 definiremos dos vectores q1 y q2, y luego los multiplicaremos
q1=[1 2 3] q2=[-1 0 0 1] conv(q1,q2)
q1 =
1 2 3
q2 =
-1 0 0 1
ans =
-1 -2 -3 1 2 3
Ej 1. Halle el producto de los dos polinomios p1(x)=x^2 + 2*x + 1, p2(x)= x^3 + x^2 + 1. Evalue el resultado en x=2. (Sol. p1(x)*p2(x)= x^5 + 3*x^4 + 3*x^3 + 2*x^2 + 2*x + 1. p1(2)*p2(2) = 117)
Polinomio de interpolación de Lagrange
Dados los nodos x0, x1, x2, x3, ..., xN y los valores en dichos nodos y0, y1, y2, y3, ..., yN , queremos hallar el polinomio
p(x)=aN*x^N + ... + a2*x^2 + a1*x +a0
tal que p(x0)=y0, p(x1)=y1, ..., p(xN)=yN.
Decir que el polinomio p(x)=aN*x^N + ... + a2*x^2 + a1*x +a0 verifica las condiciones de interpolación p(x0)=y0, p(x1)=y1, ..., p(xN)=yN equivale a decir que los coeficientes aN,...,a2,a1,a0 del polinomio p(x) satisfacen un sistema de ecuaciones lineales cuyo vector de términos independientes es el vector y=[y0; y1; y2; y3; ...; yN], y cuya matriz de coeficientes es la matriz de Vandermonde formada por los elementos x(i,j)=xi^(n-j) con i,j=0,...,n (que puede obtenerse en Matlab mediante el comando vander(x), donde x=[x0, x1, x2, x3, ..., xN]).
Por tanto, los coeficientes del polinomio de interpolación de Lagrange p(x) se pueden obtener resolviendo dicho sistema de ecuaciones lineales.
Como ejemplo nos planteamos encontrar el polinomio de grado menor o igual que 3 que cumple p(1)=1, p(2)=1, p(3)=2, p(4)=6. (Ref: Diez lecciones de Cálculo Numérico. J. M. Sanz Serna.)
x=[1, 2, 3, 4] y=[1; 1; 2; 6] vander(x)
x =
1 2 3 4
y =
1
1
2
6
ans =
1 1 1 1
8 4 2 1
27 9 3 1
64 16 4 1
Obtenemos la solución del sistema vander(x)*p=y mediante la 'división izquierda'. Para ello el vector y debe ser un vector columna.
p=vander(x)\y
p =
0.3333
-1.5000
2.1667
0.0000
Representamos el polinomio p evaluándolo en los valores del vector xd
xd=[0:0.1:5]; yd=polyval(p,xd); plot(xd,yd)
Para observar gráficamente el polinomio de interpolación representaremos los datos de interpolación (los vectores x e y) mediante círculos, y representaremos el polinomio p mediante los vectores xd e yd.
plot(x,y,'o',xd,yd)
Ej 2. Halle el polinomio q de grado cuatro que verifica que q(2)= 3, q(4)=5, q(5)=8, q(7)=10, q(10)= 4. Halle el determinante y el número de condición (usando la norma 2) de la matriz de Vandermonde asociada a los valores dados. (Observe que el determinante es distinto de cero y que el número de condición es mucho mayor que uno) (Sol. Coeficientes de q = [0.0347; -0.8917; 7.5931; -23.7583; 26.7222]. Determinante = 129600. Número de condición = 4.1024e+005)
Diferencias divididas
A continuación vamos a construir con Matlab la tabla de diferencias divididas para el polinomio de interpolación de Lagrange p(x) con nodos x(0), x(1), ..., x(N-1), x(N) (Ref: Diez lecciones de Cálculo Numérico. J. M. Sanz Serna.) (Ref: Numerical methods using Matlab. J. H. Mathews y K. D. Fink.)
Construiremos una tabla como la que se muestra, donde xN-1 es el nodo x(N-1), xN-2 el nodo x(N-2), etc.
x0 p[x0] x1 p[x1] p[x0,x1] x2 p[x2] p[x1,x2] p[x0,x1,x2] . . . . . . . . . . . . . . . xN p[xN] p[xN-1,xN] p[xN-2,xN-1,xN] ... p[x0,...,xN]
Definimos los nodos, el número de nodos N y los valores del polinomio p(x) en tales nodos.
nodos=[1 2 3 4]; N=length(nodos); pnodos=[1, 1, 2, 6];
Inicializamos una matriz M, con NaN, donde iremos almacenando las diferencias divididas.
M=NaN(N,N+1) % Rellenamos la primera columna M(:,1)=nodos % Rellenamos la segunda columna (son las diferencias divididas de orden 0) M(:,2)=pnodos
M =
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
NaN NaN NaN NaN NaN
M =
1 NaN NaN NaN NaN
2 NaN NaN NaN NaN
3 NaN NaN NaN NaN
4 NaN NaN NaN NaN
M =
1 1 NaN NaN NaN
2 1 NaN NaN NaN
3 2 NaN NaN NaN
4 6 NaN NaN NaN
Calcularemos las diferencias divididas de orden i=1, 2, ..., N-1 que colocaremos en la columna i+2 de la matriz.
for j=3:N+1 for i=j-1:N % Las diferencias divididas de orden i (situadas en la columna i+2) % comienzan en la fila i+1 y acaban en la fila N. % % Reproducimos la parte de la matriz M necesaria para hallar el elemento % M(i,j) % % M(i-j+2,1) M(i-j+2,2) % . % . % . % M(i-1,j-1) % M(i,1) . . . . . M(i,j-1) M(i,j) % % Puede comprobarse que siguiendo el elemento M(i-j+2,2) en diagonal (esto % es, sumándole j-2 al indice de las filas y de las columnas) se llega al % elemento M(i,j) . M(i,j)=(M(i,j-1) - M(i-1,j-1)) / (M(i,1) - M(i-j+2,1) ); end end
La matriz M contiene la tabla de diferencias divididas (Observemos que las diferencias divididas que se utilizan para escribir la forma de Newton del polinomio de Lagrange p(x) son los coeficientes M(i,i+1) de la matriz M con i=1,..N).
M
M =
1.0000 1.0000 NaN NaN NaN
2.0000 1.0000 0 NaN NaN
3.0000 2.0000 1.0000 0.5000 NaN
4.0000 6.0000 4.0000 1.5000 0.3333
Ej 3. Construya la tabla de diferencias divididas correspondiente al polinomio p(x) que interpola la función coseno en los nodos 0, 1 ,2, 3, 4. Localice el valor de p[0, 1 ,2] en la matriz construida. (Sol. Diferencias divididas = 1.0000 , -0.4597 , -0.2484 , 0.1466 , -0.0147) (Ref: Numerical methods using Matlab. J. H. Mathews y K. D. Fink.)
Forma de Newton del polinomio de interpolación de Lagrange
Con las diferencias divididas calculadas podemos construir el polinomio de interpolación usando la forma de Newton, donde xN-1 es el nodo x(N-1): P(x) =p[x0] + p[x0,x1]*(x-x0) + p[x0,x1,x2]*(x-x0)*(x-x1) + ... + p[x0,x1,...,xN]*(x-x0)*(x-x1)* ... *(x-xN-1) .
En primer lugar veamos cómo construir el polinomio q(x)=(x-x0)*(x-x1)* ... *(x-xN) a partir de los valores x=[x0, x1, x2, ..., xN]. Para ello, usaremos un acumulador de producto, que comenzará valiendo el polinomio constante [1], y le iremos multiplicando los monomios (x-xj) que se expresan en forma de vectores como [1, -xj].
q=[1] for i=2:N monomio=[1, -nodos(i-1)] q=conv(q,monomio) end
q =
1
monomio =
1 -1
q =
1 -1
monomio =
1 -2
q =
1 -3 2
monomio =
1 -3
q =
1 -6 11 -6
El polinomio p(x) se obtiene, usando la forma de Newton, mediante una suma de polinomios cuyos sumandos son los polinomios q (calculados como en el bucle anterior) multiplicados por la correspondiente diferencia dividida. Observe en el siguiente bucle cómo se calcula el polinomio p(x) mediante el proceso descrito utilizando un acumulador de suma:
q=[1] newton=[pnodos(1)] for i=2:N monomio=[1, -nodos(i-1)] q=conv(q,monomio) newton=[0, newton] + M(i,i+1)*q end
q =
1
newton =
1
monomio =
1 -1
q =
1 -1
newton =
0 1
monomio =
1 -2
q =
1 -3 2
newton =
0.5000 -1.5000 2.0000
monomio =
1 -3
q =
1 -6 11 -6
newton =
0.3333 -1.5000 2.1667 0
Ej 4. Escribe el polinomio p(x) que se acaba de obtener. ¿Respecto de qué base está expresado? Expresa p(x) respecto de la base de Newton del problema.
Ej 5. ¿Por qué se necesita un 0 en el sumando [0, newton]? (Observe los grados de los polinomios que se suman en la forma de Newton)
Ej 6. Si nodos es [x0, x1, ..., xN] ¿para i=2,..,N qué nodos son nodos(i-1) y qué diferencias divididas son M(i,i+1) ? (Piense el caso i=2, luego el caso i=3, ...)
Ej 7. Se quiere aproximar la función seno en el intervalo [0 , 3.14] mediante un polinomio de grado 6. i) Halle la tabla de diferencia divididas correspondiente a los nodos 0, 0.5, 1, 1.5, 2, 2.5, 3 y a los valores de la función seno. (Sol. Diferencias divididas = 0 , 0.9589 , -0.2348 , -0.1182 , 0.0336, 0.0025, -0.0013 ) ii) Halle el polinomio de interpolación correspondiente a los datos descritos usando la forma de Newton. iii) Represente conjuntamente los datos de interpolación (en rojo) y el polinomio obtenido (en azul). iv) Represente conjuntamente el polinomio obtenido (en azul) y la función seno (en rojo).
% Matemáticas II. Grado en Ingeniería Química. (A. Palomares)
disp(date)
22-Mar-2017