function A6_iii

n = 20;
x = ( cos(((2*(0:n)+1)/(2*n+2)*pi)) + 1)/2; % Tschebyscheff-Stützstellen
y = exp(x);                                 % Funktionswerte
y_eps = y + randn(size(x))/1000;            % Gestörte Werte

N  = 200;
xx = linspace(0,1,N);                       % Auswertungsstellen
yy = exp(xx);

plot(xx,yy,'k');                            % Originale e-Funktion plotten
hold on;

plot(x,y,'.k');                             % Ungestörte Werte plotten
plot(x,y_eps,'.r');                         % Gestörte Werte plotten
legend('f(x)=e^x','y_i','y_i^\epsilon',4);
drawnow;

fprintf('\nAuswertungsstellen: %i\n',N);
fprintf('Stützstellen:       %i\n\n',length(x));
fprintf('Bitte warten...\n\n');

tic; flops(0);
P = Neville_Aitken(x,y,xx);                 % Ungestörte Interpolation
fprintf('Ungestörte Stützstellen:\n');
fprintf('Zeit: %0.2fs    Flops: %i\n\n',toc, flops);
plot(xx,P,'g');
legend('f(x)=e^x','y_i','y_i^\epsilon','P mit P(x_i)=y_i',4);
drawnow;

tic; flops(0);
P_eps = Neville_Aitken(x,y_eps,xx);         % Gestörte Interpolation
fprintf('Gestörte Stützstellen:\n');
fprintf('Zeit: %0.2fs    Flops: %i\n\n',toc, flops);
plot(xx,P_eps,'r');
legend('f(x)=e^x','y_i','y_i^\epsilon','P mit P(x_i)=y_i','P^\epsilon mit P^\epsilon(x_i)=y_i^\epsilon',4);

title('Polynominterpolation mit gestörten und ungestörten Funktionswerten');
xlabel x;
ylabel('e^x');

V = max(abs(P-P_eps)) / max(abs(y-y_eps))
lambda_20 = aufgabe6_i(1,20,0)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function f = Neville_Aitken(x,y,xj)
% Neville-Aitken-Algorithmus für die Stützstellen x
% und den Funktionswerten y. Dieser Algorithmus wertet
% das Interpolationspolynom durch die Punkte (x,y)
% an den Stellen xj aus.
%
% Eingabe: x  Stützstellen
%          y  Funktionswerte, die interpoliert werden
%          xj Stellen, an denen das Polynom ausgewertet wird
%
% Ausgabe: f  = P(xj), interpolierte Werte

n = length(x);
f = zeros(size(xj));

for j=1:length(xj)
   P = y;
   
   for k=2:n
      for i=n:-1:k
         P(i) = P(i) + (xj(j)-x(i)) / (x(i)-x(i-k+1)) * (P(i)-P(i-1));
      end;
   end;
   
   f(j) = P(n);
end;

