function P = aufg01f05(matnr)

%AUFG01F05 berechnet nahezu Matrikelnummern interpolierende Polynome.
%
%   P = AUFG01F05(MATNR) berechnet mittels der Newtonschen Form
%   der Polynominterpolation die eindeutigen Interpolationspolynom
%   zu den Aufgaben: Berechne p, q so dass
%
%              p(i) = i-te Ziffer der Matrikelnummer
%              q(i) = p(i), i != 3, q(3) = p(3)+1/100.
%
%   gilt. Dabei ist die Laenge der Matrikelnummer nicht a priori
%   beschraenkt. Implizit wird die Minimallaenge 5 verwendet. Das
%   bedeutet, dass eine kuerzere Matrikelnummer fuehrende Nullen
%   besitzt, also MATNR = 911 => MATNR = 00911 .
%
%   Der ausgegebene Array P enthaelt die Vektoren der Koeffizienten
%   der beiden Polynome in der Form [P,Q].
%
%   AUFG01F05 ist die Loesung der Aufgabe 3c) von Uebungsblatt 1
%   der Vorlesung 'Numerische Verfahren' von Jens Zemke. Die Funktion
%   arbeitet unter Matlab und Octave.
%
%   SEE ALSO: aufg01f03, aufg01f04, polyval, polyfit

%   $Author: zemke $
%   $Revision: 1.4 $
%   $Date: 2004/04/15 10:40:36 $

if matnr < 0
   error('Es wurde eine negative Matrikelnummer eingegeben');
end

if matnr ~= floor(matnr)
   error('Es wurde eine nicht ganzzahlige Matrikelnummer eingegeben');
end

%
% Berechne die Laenge der Matrikelnummer, i.e. die Anzahl der Ziffern
% und transformiere die ganze Zahl in eine Zifferndarstellung.
%

n = max(floor(log10(matnr))+1,5);
x = 1:n; y = zeros(1,n);
potenz = 10.^(0:n-1);

for i = 1:n
   y(i) = floor(matnr/potenz(n-i+1));
   matnr = matnr - y(i)*potenz(n-i+1);
end

%
% Verdoppeln der Daten zu den einzelnen Punkten, setze dabei
% den dritten Wert der zweiten Zeile auf den alten Wert plus
% ein Hundertstel ... Damit kann man dann eine Schleife
% statt zweier verwenden ;-)
%

y = [y;y]; y(2,3) = y(2,3)+1/100; y = y';

%
% Berechnung der dividierten Differenzen gemaess Algorithmus 2.14
% des Skriptes 'Numerische Verfahren' von Jens Zemke.
%

h = zeros(n,2);
p = zeros(n,2);

for k = 1:n
   h(k,:) = y(k,:);
   if k > 1
      for i = k-1:-1:1
         h(i,:) = (h(i+1,:)-h(i,:))./(x(k)-x(i));
      end
   end
   p(k,:) = h(1,:);
end

%
% Generieren eines Vektors von Werten in [0,50] und berechnen der
% zugehoerigen Funktionswerte, berechnen des Fehlers, anschliessend
% plotten.
%

t = linspace(0,50,100);
tt = [t;t]; tt = tt';

pt = ones(100,1)*p(n,:);

for i = n-1:-1:1
   pt = pt.*(tt-x(i))+ones(100,1)*p(i,:);
end

f = pt(:,2)-pt(:,1);

plot(t,f);
grid('on');
title(['Das Fehlerpolynom ausgewertet auf dem Intervall [1,50]']);

%
% Wen's interessiert: der relative Fehler ;-)
%
% semilogy(t,abs(f./pt(:,1))); grid('on');
%
