function A13_1

%Anzahl Intervalle
n=8;
%Daten einlesen
load g36.dat;
%Datenbereich in n Intervalle aufspalten
[N,X]=hist(g36,n);
% N enthält die absoluten Häufigkeiten
% X den Mittelpunkt des Intervalls mit der Breite delta
delta=X(2)-X(1);
%relative Häufigkeit pro Intervallbreite
Y=N/sum(N)/delta;
%als Histogramm plotten
figure;
hold on;
bar(X,Y,1,'y');
% so, jetzt geht's weiter mit der Splineberechnung ...

r = 3;
% X (Intervallmittelwerte) auf x (Randpunkte der Intervalle) umrechnen
x = [X-delta/2 X(end)+delta/2];

% 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;
   
% (H) - Randbedingungen in Matrix einfügen
for j=0:r   M(2+n,  j+1) = spline_p(j,  x,x(  1),1);
            M(3+n,  j+1) = spline_p(j,  x,x(end),1); end;
for j=1:n-1 M(2+n,r+j+1) = spline_q(j,r,x,x(  1),1);
            M(3+n,r+j+1) = spline_q(j,r,x,x(end),1); end;
   
% Rechte-Hand-Seite aufbauen
RHS = [0;cumsum(Y')*delta;0;0];   % Interpolationsbedingungen
      
% 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
xx = linspace(min(x),max(x),1000);
S = ComputeSpline(n,r,x,alpha,beta,xx,1);

hold on;
plot(xx,S,'r','linewidth',2);
plot(x,0,'*k');



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