% CVICENI - 8
% Interpolace polynomem a splinem.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1) Interpolace algebraickym polynomem
% Tabulka dat [x,y]:
x = 1:1:10
y = [-2,-5,-3,1,4,6,8,6,7,6]
% p(x) = a0 + a1*x + a2*x^2 + ...
% ziskame resenim rovnice: A*a = y
% A(i,j) = x(i)^(j-1) ... Srovnej s metodou nejmensich ctvercu!

for i=1:length(x)
    for j=1:length(x)
        A(i,j) = x(i)^(j-1);
    end
end

% koeficienty polynomu
a = A\y'

% vyhodnoceni polynomu
xp = linspace(0,11,1000);
yp = polyval(flip(a),xp);

fig = figure(1);
hold on;
plot(x,y,'bo','MarkerFaceColor', 'b');
plot(xp,yp,'-k');
grid on;
title('Interpolace algebraickym polynomem')
ylim([min(y)-3, max(y)+8])
hold off;

% 2) Interpolace Lagrangeovym polynomem

x = 1:1:10;
y = [-4,-5,-3,5,7,6,8,2,7,6];

N = 1000; % jemnost vyjadreni interpolacniho polynomu
xp = linspace(0,11,N);
% Ln = y0*L0 + y1*L1 + ... + yn*Ln

% tvorba Lagrangeova polynomu
for i=1:length(x)
    L(i,:) = ones(1,N);
    for j=1:length(x)
        if (i ~= j)            
            L(i,:) = L(i,:) .* ((xp-x(j)) / (x(i)-x(j)));
        end
    end
end
  
% vypocet finalniho interpolacniho polynomu 
% (lze vclenit do predchozich cyklu)
yp = zeros(1,N);
for i=1:length(x)
    yp = yp + y(i) * L(i,:);
end

fig = figure(2);
hold on;
plot(x,y,'ro','MarkerFaceColor', 'r');
plot(xp,yp,'-k');
grid on;
title('Interpolace Lagrangeovymi polynomy')
ylim([min(y)-3, max(y)+15])
hold off;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3) Interpolace Nevilleovym algoritmem

x = [1, 2, 3, 4];
y = [5, 3, 8, 4];

% Interpolace v bode napr. X = 2.5
interpolated_value = NEV(2.5, x, y);
fprintf('Interpolovana hodnota v x = 2.5 je: %.4f\n', interpolated_value);

xp = linspace(0,5,1000);
for i=1:length(xp)
    yp(i) = NEV(xp(i), x, y);
end

fig = figure(3);
hold on;
plot(x,y,'go','MarkerFaceColor', 'g');
plot(xp,yp,'-k');
grid on;
title('Interpolace, Nevilluv algoritmus')
ylim([min(y)-15, max(y)+15])
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4) Interpolace Nevilleovym algoritmem




















%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ukazka rekurze
function P = NEV(X, x, y)
    n = length(x);    
    if n == 1
        P = y(1);
    else
        P1 = (X - x(1)) * NEV(X, x(2:end),   y(2:end)  );
        P2 = (x(n) - X) * NEV(X, x(1:end-1), y(1:end-1));
        P  = (P1+P2) / (x(n) - x(1));
    end
end
