%% 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.

%% 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]
%%
% 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)

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

%% 
% 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]

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

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

%%
% 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

%%
% 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

%%
% 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

%%
% 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

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