%% Blatt 11 - Programmieraufgaben:
%% Tschebyscheff-Interpolation (4 Punkte)
% Vorname Nachname Mat.-Nr.
% Gegeben seien m Tschebyscheff-Knoten
% x_i = \cos\Big(\frac{i-1}{m-1} \pi\Big)  für  i = 1, 2, \ldots, m,
% und m zugehörige Knotenwerte y_1, y_2, \ldots, y_m \in \mathbb{R}. Wie wir in den Programmieraufgaben von Blatt 5 gesehen haben, entstehen für die Berechnung des Interpolationspolynoms
% p(x) = \sum_{j=1}^{m} y_j L_j(x) .
% bei direkter Umsetzung der Lagrange-Polynome mit der Formel
% L_j(x) = \prod_{i=1 . i \neq j}^m \,\, \frac{x - x_i}{x_j - x_i}
% bei grossen Werten von m Probleme: Wegen den zu dichten Knoten x_j und der vielen Faktoren wird der Hauptnenner bei dieser Auswertung sehr klein. Unabhängig davon hat dieses direkte Vorgehen für eine Lagrange-Polynominterpolation-Auswertung bei n Punkten einen Gesamtaufwand von \mathcal{O}(n m^2) und ist damit auch deshalb äussert unvorteilhaft.
% Für allgemeine Knotenwahlstrategien (mit quadratischer Randknotendichte!) greift man grundsätzlich auf die Newton-Polynominterpolation mit der Horner-Schema-Auswertung zurück, wobei der Gesamtaufwand \mathcal{O}(m^2 + n m) beträgt. Dies ist für die Tschebyscheff-Interpolation allerdings noch immer suboptimal.
% Die Auswertung der durch die Tschebyscheff-Knoten definierten Lagrange-Polynome kann mittels einer algorithmisch sehr vorteilhaften Formel realisiert werden, die von Salzer (1972) stammt. Es lässt sich nämlich zeigen, die Lagrange-Polynome durch die Tschebyscheff-Knoten der Formel
% L_j(x) = \frac{\sum_{i=1}^{m} \frac{a_i \, \delta_{i, j}}{x - x_i}}{\sum_{i=1}^{m} \frac{a_i}{x - x_i}}   für  x \notin \{x_1, x_2, \ldots, x_m\}  mit  a_i = \begin{cases} (-1)^{i-1} & \text{für } i = 2, \ldots,  m-1 \\ \frac{(-1)^{i-1}}{2} & \text{für } i = 1 \text{ und } i = m \end{cases}   und  \delta_{i, j} = \begin{cases} 1 & \text{für } i = j \\ 0 & \text{für } i \neq j \end{cases}
% genügen. Es ergibt sich also
% p(x) = \frac{\sum_{i=1}^{m} \frac{a_i \, y_i}{x - x_i}}{\sum_{i=1}^{m} \frac{a_i}{x - x_i}}   für  x \notin \{x_1, x_2, \ldots, x_m\}.
% Diese Formel wird die baryzentrische Formel genannt. Offensichtlich hat eine Auswertung dieser Formel für x \notin \{x_1, x_2, \ldots, x_m\} einen Aufwand \mathcal{O}(m) und ist damit äusserst schnell; andernfalls ist x = x_j, was ebenfalls mit \mathcal{O}(m) Aufwand geprüft werden kann, und man setzt direkt p(x) = y_j. Eine Auswertung bei n Punkten hat damit einen optimalen Gesamtaufwand von \mathcal{O}(n m).
% Schreiben Sie hier am Ende dieses Live Scripts eine MATLAB-Funktion
% ws = lagrangeCheby(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 Auswertungspunkte 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 }  mit  x_i = \cos\Big(\frac{i-1}{m-1} \pi\Big),   \mathbf{y} = \pmatrix{ y_1 \cr y_2 \cr \vdots \cr y_m }  und  \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_m) } ,
% der mit der baryzentrischen Formel berechnet wird, zurückgibt.
%
% Mit dieser Funktion behandeln wir nun das nachfolgende Beispiel und betrachten die als function handle f definierte, k-1 mal stetig diffenenzierbare Funktion für k=6. Für m = 5, 15, \ldots, 55 berechnen und plotten Sie für jedes m in einem Subplot die echte und die Tschebyscheff-interpolierte Funktion übereinander. Vervollständigen Sie dazu den untenstehenden Codeblock. Eine in den Grundzügen sinnvolle Interpolation ergibt sich erst für m = 35, aber schon für m = 55 ist die Interpolation optisch kaum von der Funktion zu unterscheiden:

k = 6;
f = @(x) abs(exp(.1*x) .* cos(3*pi*x)) .* (exp(.1*x) .* cos(3*pi*x)).^(k-1);
zs = linspace(-1, 1, 2^11+1).';  %% Auswertungspunkte
figure();                        %% neues Bild

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

    xs = % TODO: berechnen
    ys = % TODO: berechnen
    ws = lagrangeCheby(xs, ys, zs);

    plot(xs, ys, 'ko', zs, f(zs), 'r-.', zs, ws, 'b-', 'MarkerSize', 2);
    title("m = "+num2str(m))
end

% Aufgrund der Tatsache, dass die k-1-te Ableitung Lipschitz-stetig ist, sollte sich im folgenden Konvergenztest die Rate m^{-k} erzeigen:

z = 1/7;
fz = f(z);

ms = (1:27:7*666).';
es = zeros(size(ms));

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

figure();
subplot(1, 2, 1);
semilogy(ms, es, 'b-', ms, isnan(es), 'cx', ms, ms.^(-k), 'r-');
legend('Fehler', 'NaN', ['Rate ' num2str(k)], 'Location', 'southwest');
subplot(1, 2, 2);
loglog(ms, es, 'b-', ms, isnan(es), 'cx', ms, ms.^(-k), 'r-');
legend('Fehler', 'NaN', ['Rate ' num2str(k)], 'Location', 'southwest');

%
% Funktionen:

% TODO: implementieren