Práctica 5. Interpolación
En esta práctica estudiaremos cómo definir y utilizar polinomios en Matlab o Octave, 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/Octave
En Matlab y en Octave 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 con datos 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 f=[f0, f1, f2, f3, ..., fN] buscamos un polinomio p(x)=aN*x^N + ... + a2*x^2 + a1*x +a0 de forma que p(x0)=f0, p(x1)=f1, ..., p(xN)=fN.
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 f, 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 Li(pj)=pj(xi)=xi^j. En Matlab el vector p contiene los coeficientes ordenados de mayor grado a menor grado, y la matriz del sistema, adecuada a ese orden, se puede obtener mediante el comando vander(x).
Como ejemplo nos planteamos encontrar el polinomio de grado 3 que cumple p(-1)=-16, p(0)=-7, p(1)=4, p(2)=-1.
x=[-1, 0, 1, 2] f=[-16; -7; -4; -1] vander(x)
x =
-1 0 1 2
f =
-16
-7
-4
-1
ans =
-1 1 -1 1
0 0 0 1
1 1 1 1
8 4 2 1
Obtenemos la solución del sistema vander(x)*p=f mediante la 'división izquierda'. Para hacer esta división el vector f debe ser un vector columna.
p=vander(x)\f
p =
1
-3
5
-7
Representamos el polinomio p evaluándolo en los valores del vector xd .
xd=[-2:0.1:3]; 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 y f) mediante círculos, y representaremos el polinomio p mediante los vectores xd e yd.
plot(x,f,'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. (Sol. Coeficientes de q = [0.034722; -0.891667; 7.593056; -23.7583; 26.7222]. )
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).
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=[-3 -2 -1 0 2]; N=length(nodos)-1; fnodos=[23, -3, -3, -1, 33];
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 =
-3 23 NaN NaN NaN NaN
-2 -3 -26 NaN NaN NaN
-1 -3 0 13 NaN NaN
0 -1 2 1 -4 NaN
2 33 17 5 1 1
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 cómo construir el polinomio q(x)=(x-x0)*(x-x1)* ... *(x-xN-1) a partir de los valores x=[x0, x1, x2, ..., xN-1]. 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 3
q =
1 3
monomio =
1 2
q =
1 5 6
monomio =
1 1
q =
1 6 11 6
monomio =
1 0
q =
1 6 11 6 0
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 3
q =
1 3
newton =
-26 -55
monomio =
1 2
q =
1 5 6
newton =
13 39 23
monomio =
1 1
q =
1 6 11 6
newton =
-4 -11 -5 -1
monomio =
1 0
q =
1 6 11 6 0
newton =
1 2 0 1 -1
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 el vector [x0, x1, ..., xN] ¿qué nodo entre x0, x1, ..., xN 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.
% Ecuaciones Diferenciales y Cálculo Numérico. (A. Palomares) % Grado en Ingeniería de Tecnologías de Telecomunicación. disp(date)
24-Apr-2023