%% Blatt 5 - Programmieraufgaben:
%% Lagrange-Interpolation (4 Punkte)
% Vorname Nachname Mat.-Nr.
% Gegeben seien m verschieden Knoten x_1, x_2, \ldots, x_m \in \mathbb{R} und m zugehörige Knotenwerte y_1, y_2, \ldots, y_m \in \mathbb{R}. Bekanntlich gibt es dann ein eindeutiges Polynom p \in \Pi_{m-1}, also ein Polynom vom Grad m-1, welches diese Werte interpoliert, y_i = p(x_i). Eine mögliche Darstellung dieses Polynom beruht auf den Lagrange-Polynomen, L_1, L_2, \ldots, L_m, welche durch
% L_j(x) = \prod_{i=1 . i \neq j}^m \,\, \frac{x - x_i}{x_j - x_i}
% gegeben sind. Damit ergibt sich das Interpolationspolynom als
% p(x) = \sum_{j=1}^{m} y_j L_j(x) .
%
% Schreiben Sie hier am Ende dieses Live Scripts eine MATLAB-Funktion
% ws = lagrange(xs, ys, zs) ,
% welches die Knoten x_1, x_2, \ldots, x_m \in \mathbb{R}, die Knotenwerte y_1, y_2, \ldots, y_m \in \mathbb{R} und die Auswertungsorte z_1, z_2, \ldots, z_n \in \mathbb{R} als Spaltenvektoren
% \mathbf{x} = \pmatrix{ x_1 \cr x_2 \cr \vdots \cr x_m }, \quad \mathbf{y} = \pmatrix{ y_1 \cr y_2 \cr \vdots \cr y_m } \quad\text{und}\quad \mathbf{z} = \pmatrix{ z_1 \cr z_2 \cr \vdots \cr z_n }
% als Eingabe nimmt und den Spaltenvektor
% \mathbf{w} = \pmatrix{ p(z_1) \cr p(z_2) \cr \vdots \cr p(z_n) }
% zurückgibt.
%
% Mit dieser Funktion betrachten wir nun das Beispiel von Runge: Wir betrachten für gegebene Knoten x_j die Knotenwerte y_j = f(x_j), wobei f(x) = \frac{1}{1+25x^2}. Wählen Sie die m Knoten als äquidistant von -1 bis 1, für m = 5, 9, \ldots, 25, und plotten Sie für jedes m in einem Subplot die echte und die Interpolierte Funktion übereinander:

f = @(x) 1 ./ (1 + 25*x.^2);  % Runge Funktion
zs = linspace(-1, 1, 666).';  % Auswertungspunkte
figure();                     % neues Bild

p = 0;
for m=5:4:25
    p = p+1;
    subplot(2, 3, p);         % setze aktueller Subplot

    % TODO: Implementieren
    
end

% Die Randoszillationen lassen sich durch quadratisches Anheben der Knotendichte an den Rändern, z.B. indem man die Chebyshev-Knoten,
% x_i = \cos\Big(\frac{i-1}{m-1} \pi\Big) für i = 1, 2, \ldots, m,
%  benutzt, vermeiden. Wiederholen Sie das obige Runge-Beispiel für die Chebyshev-Knoten:

figure();                     % neues Bild

p = 0;
for m=5:4:25
    p = p+1;
    subplot(2, 3, p);         % setze aktueller Subplot

    % TODO: Implementieren
    
end

% Mit den Chebyshev-Knoten ist die Interpolationsaufgabe offenbar wohlkonditioniert. Wie man anhand von nachfolgendem Beispiel sieht, erreicht die Interpolation für sehr glatte Funktion schnelle Konvergenzraten. Für k-fach stetig differenzierbare Funktionen verhält sich der Fehler beweisbar wie m^{-k}, bei analytischen Funktionen sogar exponentiell:

k = 2;
f = @(x) x.^(k-1) .* abs(x);
z = 1/7;
fz = f(z);

ms = (1:4:999).';
es = zeros(size(ms));

tic;
for j=1:length(ms)
    xs = cos(linspace(0, pi, ms(j))).';
    ys = f(xs);
    w = lagrange(xs, ys, z);
    es(j) = abs(fz - w);
end
toc

figure();
semilogy(ms, es, 'b-', ms, isnan(es), 'rx', ms, ms.^(-k-1), 'r-');
legend('Feler', 'NaN', ['Rate ' num2str(k+1)], 'Location', 'north');

% Man sieht aber, dass die direkte Auswertung der Lagrange-Polynome Probleme bereiten kann. Nach m= 800 steigt der Fehler wieder an und vor m = 1000 werden dann sogar nur noch NaN Werte ausgegeben; dies ist auf eine direkte Berechnung der Werte L_j(z) mit der Produktdarstellung zurück zu führen. In anderen Worten: Die Berechnung der L_j(z) mit Chebyshev-Knoten ist zwar wohlkonditioniert, aber die gewählte Berechnungsart instabil. Ausserdem hat diese Umsetzung den grossen Nachteil für eine Auswertung bei n Punkten einen Gesamtaufwand von \mathcal{O}(n m^2) zu besitzen. Für allgemeine Knotenwahlstrategien (mit quadratischer Randknotendichte!) sollte deshalb auf die Newton-Polynominterpolation mit der Horner-Schema-Auswertung zurückgegriffen werden, da dies dann stabil und mit Gesamtaufwand \mathcal{O}(m^2 + n m) auch effizienter ist.
%
% Funktionen:

function ws = lagrange(xs, ys, zs)
    % TODO: Implementieren
    
end