function A12(Randbedingung, Funktion)
% A12(Randbedingung, Funktion)
%
% Untersuchung des Konvergenzverhaltens von Splines
% Eingabe: Randbedingung:
%              0 für Hermite-Splines
%              1 für natürliche Splines
%          Funktion:
%              0 für cos(2*pi*x)
%              1 für |x|^3

n  = 6;
m  = 1;

x  = linspace(0,1,n+1);  % Stützstellen
xx = linspace(0,1,500);  % Auswertungsstellen
r  = 2*m+1;              % Grad der Splines

if 1                     % Zum Zeichnen der Basis-Splines auf "1" setzen
   figure;
   plot(x,0,'*r');
   hold on;
   for j=0:r
      plot(xx,spline_p(j,x,xx,0),'k');
   end;
   
   for j=1:n-1
      plot(xx,spline_q(j,r,x,xx,0),'r');
   end;
   
   axis([x(1) x(end) -0.2 1.2]);
   xlabel x
   title('Basissplines p_j (schwarz) und q_j (rot)');
   
   figure;
   plot(xx,spline_q(2,r,x,xx,0),'k');
   hold on;
   plot(xx,spline_q(2,r,x,xx,1),'--k');
   plot(xx,spline_q(2,r,x,xx,2),'-.k');
   plot(xx,spline_q(2,r,x,xx,3),':k');   
   plot(x,0,'*r');
   
   xlabel x
   title('Ableitungen der Basissplines q_2');
   legend('q_2','d/dx q_2','d^2/dx^2 q_2','d^3/dx^3 q_2');
   axis([x(1) x(end) -0.2 6.2]);   
end;

% Matrix aufsetzen

flops(0); tic;

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);

fprintf('\nAufstellen der Matrix, Berechnung der Koeffizienten\n');
fprintf('und Auswerten des Splines (%i Stellen):\n',length(xx));
fprintf('Zeit = %0.2fs, Flops: %i\n', toc, flops);

figure;
spy(M);

figure;
plot(xx,f(xx,0,Funktion),'k');
hold on;
plot( x,f( x,0,Funktion),'*r');
plot(xx,S,'r');
xlabel x
title('Splineinterpolation');
legend('Funktion f','Stützstellen','Spline');

% Numerische Integrationsstellen setzen
xxx    = (0.5:16383.5)/16384;
h      = xxx(2)-xxx(1);

% Biegeenergie und linearisierte Biegeenergie numerisch integrieren
fE_LIN = sum(f(xxx,2,Funktion).^2)*h;
fE     = sum(f(xxx,2,Funktion).^2 ./ (1+f(xxx,1,Funktion).^2).^3 )*h;

sE_LIN = sum(ComputeSpline(n,r,x,alpha,beta,xxx,2).^2)*h;
sE     = sum(ComputeSpline(n,r,x,alpha,beta,xxx,2).^2 ./ (1+ComputeSpline(n,r,x,alpha,beta,xxx,1).^2).^3 )*h;

fprintf('\nLinearisierte Biegeenergie von f=%11.6f, von S=%11.6f',fE_LIN,sE_LIN);
fprintf('\n              Biegeenergie von f=%11.6f, von S=%11.6f\n\n',fE,sE);



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,func)
%
% Wertet die order-ste Ableitung der Funktion f
% an den Stellen x aus.
%
if func==0
   Faktor = 2*pi;
   if order==0 ff =  cos(xx*Faktor);
   elseif order==1 ff = -sin(xx*Faktor);
   elseif order==2 ff = -cos(xx*Faktor);
   elseif order==3 ff =  sin(xx*Faktor);
   elseif order==4 ff =  cos(xx*Faktor);
   elseif order==5 ff = -sin(xx*Faktor);
   elseif order==6 ff = -cos(xx*Faktor);
   elseif order==7 ff =  sin(xx*Faktor);
   else error('NOT IMPLEMENTED!!!'); end;

   ff = ff * Faktor^order;
elseif func==1
   ff = spline_q(0,8,0.5,xx,order);
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;   
