function A15(Randbedingung, Funktion)
% A15(Randbedingung, Funktion)
%
% Untersuchung des Konvergenzverhaltens von Splines
% Eingabe: Randbedingung:
%              0 für Hermite-Splines
%              1 für natürliche Splines
%          Funktion:
%              0 für 1/(1+25x^2)
%              1 für |x|^3

xx   = linspace(-0.5,0.5,500);  % Auswertungsstellen
nVec = 11:2:59;
hVec = 1./nVec;
err  = hVec;
flops(0); tic;

for Count=1:length(nVec)
   n = nVec(Count);   
   m = 1;
   
   fprintf('n = %3i\n',n);
   
   x = linspace(-0.5,0.5,n+1);  % Stützstellen
   r = 2*m+1;                   % Grad der Splines
   
   % Matrix aufsetzen   
   
   M = zeros(n+r);
   for j=0:r
      M(1:n+1,j+1) = spline_p(j,x,x,0)';
   end;
   
   for j=1:n-1
      M(1:n+1,r+j+1) = spline_q(j,r,x,x,0)';
   end;
   
   if Randbedingung==0
      % (H) - Randbedingungen in Matrix einfügen
      for order=1:m
         for j=0:r
            M(order+n+1  ,j+1) = spline_p(j,x,x(  1),order);
            M(order+n+1+m,j+1) = spline_p(j,x,x(end),order);
         end;
         
         for j=1:n-1
            M(order+n+1  ,r+j+1) = spline_q(j,r,x,x(  1),order);
            M(order+n+1+m,r+j+1) = spline_q(j,r,x,x(end),order);
         end;
      end;
   else
      % (N) - Randbedingungen in Matrix einfügen
      for order=m+1:2*m
         for j=0:r
            M(order+n+1-m,j+1) = spline_p(j,x,x(  1),order);
            M(order+n+1  ,j+1) = spline_p(j,x,x(end),order);
         end;
         
         for j=1:n-1
            M(order+n+1-m,r+j+1) = spline_q(j,r,x,x(  1),order);
            M(order+n+1  ,r+j+1) = spline_q(j,r,x,x(end),order);
         end;
      end;
   end;
   
   % Rechte-Hand-Seite aufbauen
   RHS = zeros(n+r,1);
   RHS(1:n+1) = f(x,0,Funktion)';   % Interpolationsbedingungen
   
   if Randbedingung==0
      % (H) - Rechte-Hand-Seite einfügen
      for order=1:m
         RHS(order+n+1  ) = f(x(  1),order,Funktion);
         RHS(order+n+1+m) = f(x(end),order,Funktion);
      end;
   end;
   % (N) - Rechte-Hand-Seite muss nicht explizit
   % bestimmt werden, da sie sowieso Null ist.
   
   % Koeffizienten berechnen
   Coeff = M \ RHS;
   alpha = Coeff(1:r+1);      % Koeffizienten der p_j
   beta  = Coeff(r+2:end);    % Koeffizienten der q_j
   
   % Spline ausrechnen und plotten
   S = ComputeSpline(n,r,x,alpha,beta,xx,0);
         
   err(Count) = max(abs(S-f(xx,0,Funktion)));
end;

fprintf('\nZeit = %0.2fs, Flops: %i\n\n', toc, flops);

loglog(hVec,err,'.-k');
hold on;
grid on;
loglog(hVec,hVec.^1.5,'.-r');
loglog(hVec,hVec.^3,'.-g');
loglog(hVec,hVec.^4,'.-b');
xlabel h
ylabel Fehler

if (Randbedingung==0) & (Funktion==0) title('Spline-Interpolationsfehler (f_1, (H))'); end;
if (Randbedingung==1) & (Funktion==0) title('Spline-Interpolationsfehler (f_1, (N))'); end;
if (Randbedingung==0) & (Funktion==1) title('Spline-Interpolationsfehler (f_2, (H))'); end;
if (Randbedingung==1) & (Funktion==1) title('Spline-Interpolationsfehler (f_2, (N))'); end;

legend('Fehler','h^{1.5}','h^3','h^4');
axis tight;



function S = ComputeSpline(n,r,x,alpha,beta,xx,order)
S = zeros(size(xx));
for j=0:r
   S = S + alpha(j+1) * spline_p(j,x,xx,order);
end;

for j=1:n-1
   S = S + beta(j) * spline_q(j,r,x,xx,order);
end;




function ff = f(xx,order,funktion)
%
% Wertet die order-ste Ableitung der Funktion f
% an den Stellen x aus.
%
if funktion==1
   if order==0 ff =  abs(xx).^3;
   elseif order==1 ff = 3*xx.^2.*sign(xx);
   elseif order==2 ff = 6*abs(xx);
   else error('NOT IMPLEMENTED!!!'); end;
else
   if order==0 ff = 1./(1+25*xx.^2);
   elseif order==1 ff = -50*xx ./ (1+25*xx.^2).^2;
   else error('NOT IMPLEMENTED!!!'); end;
end;


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

function p = spline_p(j,x,xx,order)
% p = spline_p(j,x,xx,order)
%
% Diese Funktion wertet die order-ste Ableitung
% der Funktion p_j(x) = (x-x(1))^j an den Stellen
% xx aus. x ist der Vektor der Stützstellen.

if order>j
   p = zeros(size(xx));
else
   p = ((xx-x(1)).^(j-order));
   for n=j-order+1:j
      p=p*n;
   end;
end;

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

 
function q = spline_q(j,r,x,xx,order)
% q = spline_q(j,r,x,xx,order)
%
% Diese Funktion wertet die order-ste Ableitung der
% Funktion q_j(x) = ((x-x(j+1))^r)_+ an den Stellen
% xx aus. x ist der Vektor der Stützstellen.

q = (xx>=x(j+1)) .* (xx-x(j+1)).^(r-order) * (order<=r);
for n=r-order+1:r
   q=q*n;
end;   
