function A18;
global AuswertStellen;

AuswertStellen = 200;

tic; flops(0);
PlotteWas(1);
fprintf('\nAuswertung eines Lagrange-Splines (%i Stellen)\n',AuswertStellen);
fprintf('Zeit = %0.2fs, Flops: %i\n\n', toc, flops);

n  = 30;
r  = 3;
h  = 1/n;
x  = -r*h:h:(n+r)*h;                     % Stützstellen
xx = linspace(0,1-1e-8,AuswertStellen);  % Auswertungsstellen

tic; flops(0);
P  = psi(5,0:h:1,xx);
fprintf('Zum Vergleich: Auswertung eines Lagrange-Polynoms\n');
fprintf('Zeit = %0.2fs, Flops: %i\n\n', toc, flops);

plot(xx,P,'');
legend('Lagrange-Spline','Stützstellen','Lagrange-Polynom')

return;




for n=[0 2 3 4 5]
  PlotteWas(n);
end;


function PlotteWas(WasDenn)
global AuswertStellen;
n  = 12;
r  = 3;
h  = 1/n;
x  = -r*h:h:(n+r)*h;                     % Stützstellen
xx = linspace(0,1-1e-8,AuswertStellen);  % Auswertungsstellen

figure;
    if WasDenn==0 S = BerechneSpline(n,x,xx,(0:n)==0    ,0,0);
elseif WasDenn==1 S = BerechneSpline(n,x,xx,(0:n)==4    ,0,0);
elseif WasDenn==2 S = BerechneSpline(n,x,xx,(0:n)==8    ,0,0);
elseif WasDenn==3 S = BerechneSpline(n,x,xx,zeros(n+1,1),1,0);
elseif WasDenn==4 S = deBoor(x,r,(0:n+3)==3,xx);
elseif WasDenn==5 S = deBoor(x,r,(0:n+3)==8,xx); end;
   
plot(xx,S,'k');
hold on;
plot(x(r+1:end-r),zeros(n+1,1),'*r');
grid on;
    if WasDenn==0 title('Lagrange-Spline S_0');
elseif WasDenn==1 title('Lagrange-Spline S_4');
elseif WasDenn==2 title('Lagrange-Spline S_8');
elseif WasDenn==3 title('Lagrange-Spline T_0');
elseif WasDenn==4 title('B_{0,3}');
elseif WasDenn==5 title('B_{5,3}'); end;
drawnow;



function S = BerechneSpline(n,x,xx,f,fPrime_0,fPrime_n);
r = 3;
h = 1/n;

% Koeffizienten berechnen, d.h. erst die Matrix aufsetzen ...
A = 4*eye(n+r) + diag(ones(n+r-1,1),-1) + diag(ones(n+r-1,1),1);
A(  1,[1:r    ])=[-1 0 1];
A(n+r,[n+1:n+r])=[-1 0 1];

% ... dann die rechte Seite ...
d          = zeros(n+r,1);
d([2:n+2]) = f(:)*6;
d(1)       = 2*fPrime_0*h;
d(n+r)     = 2*fPrime_n*h;

% ... dann das Gleichungssystem lösen.
a = A \ d;

% Schliesslich noch deBoor drauf werfen und fertig ist der Spline
S = deBoor(x,r,a,xx);


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

function S = deBoor(x,r,a,xx)
% S = deBoor(x,r,a,xx)
%
% DeBoor-Algorithmus für Grad r, die Stützstellen x
% und den Koeffizientenvektor a. Dieser Algorithmus wertet
% den Spline S = sum_{i=0}^n a_i B_{i,r} an den Stellen xx aus.
%
% Eingabe: x  Stützstellen
%          r  Splinegrad  
%          a  Koeffizientenvektor (Gewichte der B_{i,r})
%          xx Stellen, an denen der Spline ausgewertet wird
%
% Ausgabe: S  = S(xx), interpolierte Werte

x = x(:);
S = zeros(size(xx));

for nStSt=1:length(xx)
   % Finde j, so dass x(j) <= xx(nStSt) <= x(j+1)
   j = min(find((xx(nStSt)>=x(1:end-1)) & (xx(nStSt)<x(2:end))));
   b = a(:);
   
   for k=1:r
      i    = j-r+k:j;
      w_ir = (xx(nStSt)-x(i)) ./ (x(i+r-k+1)-x(i));
      b(i) = w_ir .* b(i) + (1-w_ir) .* b(i-1); 
   end;
   
   S(nStSt) = b(j);
end;


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

function f = psi(i,x,xj);
% Werte die Funktion psi_i für die Stützstellen x
% an den Stellen xj aus. ACHTUNG: der erste index i
% ist 1, nicht 0 (im Gegensatz zur Vorlesung)!
%
% Eingabe: i  1 <= i <= n=length(xj)
%          x  Stützstellen
%          xj Stellen, an denen psi_i ausgewertet wird
%
% Ausgabe: f  = psi_i(xj),

f = zeros(size(xj));

for j=1:length(xj)
   f(j) = prod((xj(j)-x([1:i-1 i+1:length(x)])) ./ (x(i)-x([1:i-1 i+1:length(x)])));
end;


