%% Blatt 6 - Programmieraufgaben:
%% Newton-Interpolation (4 Punkte)
% Vorname Nachname Mat.-Nr.
% Gegeben seien m verschiedene 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 Polynoms beruht auf der Interpolation nach Newton, in welcher das Interpolationspolynom gemäss
% p(x) = \sum_{j = 0}^{n} f[x_0, \ldots, x_j] N_j(x)
% gegeben ist. Dabei sind f[x_0, \ldots, x_j] die dividierten Differenzen, welche rekursiv durch f[x_j] = y_j sowie
% f[x_j, x_{j+1}, \ldots, x_{j+i}] = \frac{f[x_{j+1}, \ldots, x_{j+i}] - f[x_j, \ldots, x_{j+i-1}]}{x_{j+i} - x_j}
% definiert sind.
% Für die Auswertung des Polynoms p sind allerdings nicht alle dividierten Differenzen notwendig, sondern nur f[x_0], \ldots, f[x_0, \ldots, x_n], welche sich leicht gemäss dem folgenden Algorthmus berechnen lassen:
% Input: Stützstellen ({\bf x, y}) \in \mathbb{R}^{n+1} \times \mathbb{R}^{n+1}
% Output: Dividierte Differenzen f[x_0], \ldots, f[x_0, \ldots, x_n], gespeichert in {\bf N} \in \mathbb{R}^{n+1}
%     setze {\bf N} := {\bf y}
%     for j = 0, \ldots, n-1 do
%     setze \mathbf{N}_{j+1 : n} := (\mathbf{N}_{j+1:n} - \mathbf{N}_{j : n-1})./(\mathbf{x}_{j+1:n} - \mathbf{x}_{0:n-j-1})
%     end for
% Beachen Sie, dass sich die Indizes bei \bf Nund \bf x auf die Indizes der Stützstellen und nicht auf den Eintrag im Vektor beziehen!
%
% Sind die dividierten Differenzen einmal berechnet, kann das Polynom p effizient mit dem Horner-Schema ausgewertet werden. Dieses lässt sich folgendermassen realisieren:
% Input: Dividierte Differenzen {\bf N} \in \mathbb{R}^{n+1}, Stützstellen {\bf x} \in \mathbb{R}^{n+1}, Auswertepunkte {\bf \xi} \in \mathbb{R}^m
% Output: Werte des Interpolationspolynoms {\bf p} \in \mathbb{R}^m
%     setze {\bf p} := {\bf N}_n*ones(size(\bf \xi) )
%     for j = n-1, \ldots, 0 do
%         setze {\bf p} := {\bf p} .*({\bf \xi} - {\bf x}_j) + {\bf N}_j
%     end for
% Dabei soll {\bf \xi} auch ein Vektor sein können.
%
% Schreiben Sie hier am Ende dieses Live Scripts MATLAB-Funktionen
% N = divDiff(x, y),
% p = hornerNewton(N, x, xi),
% welche die dividierten Differenzen berechnen, respektive, das Polynom p in Newton-Darstellung in den Punkten \bf \xi mittels Horner-Schema auswerten.
%
% Um eine glatte, parametrische Kurve durch vorgegebene Punkte {\bf p}_0 = (x_0, y_0, z_0), \ldots, {\bf p}_n = (x_n, y_n, z_n) zu legen, kann man wie folgt vorgehen: Man definiert t_0 := 0 und t_{j+1} := t_j + (x_{j+1} - x_j)^2 + (y_{j+1} - y_j)^2 + (z_{j+1} - z_j)^2, für j = 0, \ldots, n-1. Seien dann p_x, p_y und p_z die Interpolationspolynome zu den Stützstellen (t_j, x_j)_{j = 0}^n, (t_j, y_j)_{j=0}^n, respektive (t_j, z_j)_{j=0}^n, so ist
% {\bf r} : [t_0, t_n] \to \mathbb{R}^3, \qquad t \to \big( p_x(t), \, p_y(t), \, p_z(t) \big)
% eine glatte, parametrische Kurve durch die vorgegebenen Punkte. Ergänzen Sie folgenden Skelett-Code mit der Berechnung von den t_j, um eine glatte, parametrische Kurve durch die vorgegebenen Punkte zu plotten.

x = [-1.5; -1; 0; 1; 0; -1; 0; 1; 0; -1; 0; 1; 0; -1; -1.5];
y = [0; 0; -1; 0; 1; 0; -1; 0; 1; 0; -1; 0; 1; 0; 0];
z = [-6; -6; -5; -4; -3; -2; -1; 0; 1; 2; 3; 4; 5; 6; 6]./6;
% TODO: t = ...
Nx = divDiff(t, x); Ny = divDiff(t, y); Nz = divDiff(t, z);
tr = linspace(t(1), t(end), 666);
xr = hornerNewton(Nx, t, tr); yr = hornerNewton(Ny, t, tr);
zr = hornerNewton(Nz, t, tr);
plot3(xr, yr, zr, '-', x, y, z, 'ks');
axis equal; view(-35, 40); axis([-1.6, 1.1, -1.1, 1.1, -1.1, 1.1]);

%
% Funktionen:

% TODO: Implementieren