Práctica 8. Interpolación
En esta práctica estudiaremos cómo definir y utilizar polinomios en Matlab, cómo hallar el polinomio de interpolación de Lagrange, cómo calcular diferencias divididas y la forma de Newton del polinomio de interpolación.
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) podriamos hacerlo definiendo un vector de abscisas xD y evaluar 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.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 Lagrange
Dado un vector de nodos x=[x0, x1, x2, x3, ..., xN] y un vector de valores y=[y0, y1, y2, y3, ..., yN] buscamos un polinomio p(x)=aN*x^N + ... + a2*x^2 + a1*x +a0 de forma que p(x0)=y0, p(x1)=y1, ..., p(xN)=yN.
Como se vió en clase, los coeficientes del polinomio se pueden obtener resolviendo un sistema lineal de ecuaciones donde el término independiente es el vector y, el vector de incógnitas contiene los coeficientes de p(x) y la matriz del sistema es la matriz de Vandermonde formada por los elementos x(i,j)=x(i)^(n-j). La matriz del sistema se puede obtener en Matlab mediante el comando vander(x).
Como ejemplo nos planteamos encontrar el polinomio de grado 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 hacer esta división 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.034722; -0.891667; 7.593056; -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 de una función en unos 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 f[x0] x1 f[x1] f[x0,x1] x2 f[x2] f[x1,x2] f[x0,x1,x2] . . . . . . . . . . . . . . . xN f[xN] f[xN-1,xN] f[xN-2,xN-1,xN] ... f[x0,...,xN]
Definimos los nodos, N y los valores de la función.
nodos=[1 2 3 4]; N=length(nodos)-1; fnodos=[1, 1, 2, 6];
Inicializamos una matriz con NaN.
M=NaN(N+1,N+2); % Rellenamos la primera columna M(:,1)=nodos; % Rellenamos la segunda columna M(:,2)=fnodos;
Calcularemos las diferencias divididas de orden i=2, 3, ..., N+1 que colocaremos en la columna i+1 de la matriz.
for i=2:N+1 % Las diferencias divididas de orden i comienzan en la fila i y acaban % en la fila N+1. for j=i:N+1 % % Reproducimos una parte de la matriz M, necesaria para hallar el elemento % M(j,i+1) % % M(j-i+1,1) M(j-i+1,2) % . % . % . % M(j-1,i) % M(j,1) . . . . . M(j,i) M(j,i+1) % % Puede comprobarse que siguiendo el elemento M(j-i+1,2) en diagonal (esto % es, sumando i-2 al indice de las filas y de las columnas) se llega al % elemento M(j-1,i) . M(j,i+1)=(M(j,i) - M(j-1,i)) / (M(j,1) - M(j-i+1,1) ); end end
La matriz M contiene la tabla de diferencias divididas.
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 matriz de diferencias divididas correspondiente al polinomio que interpola la función coseno en los nodos 0, 1 ,2, 3, 4. Localice el valor de f[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
Con las diferencias divididas calculadas podemos construir el polinomio de interpolación en la forma de Newton, donde xN-1 es el nodo x(N-1): P(x) =f[x0] + f[x0,x1]*(x-x0) + f[x0,x1,x2]*(x-x0)*(x-x1) + ... + f[x0,x1,...,xN]*(x-x0)*(x-x1)* ... *(x-xN-1) .
En primer lugar veremos como construir el polinomio q(x)=(x-x0)*(x-x1)* ... *(x-xN) a partir de los valores x=[x0, x1, x2, ..., xN]. 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+1 monomio=[1, -nodos(i-2+1)] q=conv(q,monomio) end
monomio =
1 -1
q =
1 -1
monomio =
1 -2
q =
1 -3 2
monomio =
1 -3
q =
1 -6 11 -6
Sin embargo, la forma de Newton se expresa mediante una suma de polinomios, cuyos sumandos son los polinomios q que se calculan en el bucle, multiplicados por la correspondiente diferencia dividida. Observe en el siguiente bucle cómo la forma de Newton se calcula mediante un acumulador de suma.
newton=[fnodos(1)]; q=[1]; for i=2:N+1 monomio=[1, -nodos(i-2+1)] q=conv(q,monomio) newton=[0, newton] + M(i,i+1)*q end
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. ¿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 5. Si nodos es [x0, x1, ..., xN] ¿qué nodo es nodos(i-2+1)? ¿qué diferencia dividida es M(i,i+1) ? (Piense el caso i=2, luego el caso i=3, ...)
Ej 6. Se quiere aproximar la función seno en el intervalo [0 , 3.14] mediante un polinomio de grado 6.
i) Halle la matriz 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. ii) Halle la forma de Newton del polinomio de interpolación correspondiente a los datos descritos. iii) Represente conjuntamente los datos de interpolación y el polinomio obtenido. iv) Represente conjuntamente el polinomio obtenido y la función seno.
% Matemáticas II. Grado en Ingeniería Química. (A. Palomares)
disp(date)
28-Apr-2016
No hay comentarios:
Publicar un comentario