%Interaktiver Remez mit Einzelaustausch f\"ur f\"unf Testfunktionen:
%Fall 1: Glockenkurve exp(-x^2), x in [a,b] = [0,1],
%Fall 2: Bessel J_1, x in [0,7],
%Fall 3: xsin(1/x), x in [0.05,1]: 
%        hat 6 Nullstellen: 0.3183, 0.1592, 0.1061, 0.0796, 0.0637, 0.0531,
%Fall 4: sinus pi/4x, x in [0,1],
%Fall 5: tangens pi/4x, x in [0,1].
%Auruf: function [y,a,b]=remeztestf(j,x); [j=1,2,3,4,5] 
%Ein Schritt des Remez-Algorithmus wird ausgef\"uhrt auf die
%auf [a,b] definierte Funktion 'remeztestf' (ist separat definiert).
%Aufruf [y,a,b]=remeztestf(j,x), j ist die Funktionennummer.
%Die Knoten sind die Komponenten des Vektors Knoten: t_1,...,t_{n+1}
%Als L\"osung kommen eine obere Schranke =max |remeztestf(j,t)-Loesg(t)|, 
%t in [a,b], eine untere Schranke a0 und die
%verbesserten Knoten  heraus. Als Basis werden die Monome benutzt.
%Das resultierende Gleichungssystem
%        f(t_j)-p(t_j)=a_{n+1}(-1)^j, j=1,2,...,n+1
%mit              p(t)=a_1 + a_2t + ... + a_nt^{n-1}
%wird \"uber Interpolationstechniken (n^2 Operationen) gel\"ost.
%\"Uber Fall (1 bis 5) werden die Testfunktion und das Intervall ausgew\"ahlt.
%function [unten,oben,Knoten]=remez(Knoten,Fall);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%MATLAB-Programm, Gerhard Opfer, Approximation, SS 1995, s. remez.info
%leicht modifiziert Juni 1998, f\"ur "Numerische Mathematik f\"ur Anf\"anger"
%Wieder hervorgeholt für Approximation SS 1999
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [unten,oben,Knoten]=remez(Knoten,Fall);
n=length(Knoten)-1;         %Dimension des Raumes, Polynomgrad=n-1
z=Knoten(1,1:n);            %letzte Komponente weglassen
[f1,a,b]=remeztestf(Fall,z);%Werte an den Knoten
f2=ones(size(z));
f2(2:2:n)=-f2(2:2:n);  %Besetzung mit +- Eins
p1=polyfit(z,f1,n-1);  %Interpolationspolynom durch die ersten n Knoten
p2=polyfit(z,f2,n-1);  %dito mit Werten plus minus 1
w1=polyval(p1,Knoten(n+1));
w2=polyval(p2,Knoten(n+1));
a0=(w1-remeztestf(Fall,Knoten(n+1)))/(w2-(-1)^(n));
unten=abs(a0);         %untere Schranke f\"ur die Minimalabweichung
p=p1-a0*p2;            %Koeffizienten des gesuchten Polynoms
                       %Remezrechnung ist hier bereits zu Ende, 
xx=linspace(a,b,513);  %Zum Zeichnen
yy=remeztestf(Fall,xx);
zz=zeros(size(Knoten));
y=polyval(p,xx);
Feh=yy-y;
oben=max(abs(Feh));    %obere Schranke f\"ur die Minimalabweichung
plot(xx,Feh,[a,b],[0,0],[a,b],[a0,a0],[a,b],[-a0,-a0],Knoten,zz,'x');
t=text(-0.1,abs(a0),'');
set(t,'string','|a_0|');%Hinweis von Scott French, MATLAB
t=text(-0.11,-abs(a0),'');
set(t,'string','-|a_0|');
hold on
for j=1:n+1
  plot([Knoten(j),Knoten(j)],[0,(-1)^(j-1)*a0],'x');
  plot([Knoten(j),Knoten(j)],[0,(-1)^(j-1)*a0]);
end
Fktnamen=['exp(-x^2) ';...   %m\"ussen alle gleich lang sein!
          'Bessel J_1';...
          'xsin(1/x) ';...
          'sin \pix/4';...
          'tan \pix/4'];
xlabel(['Beste Polynomapproximation (Grad=',num2str(n-1),') von Funktion ',...
        Fktnamen(Fall,:)]);

%Wahl eines neuen Knotens:
disp('2 Punkte Anklicken: Der erste wird weggelassen, der 2. wird hinzugef\"ugt.');
disp('Ggf. wiederholen, Ende durch Enter-Taste.');
disp('Rechnung sofort beenden: x zweimal au\ss erhalb [a,b] anklicken und Enter-Taste.');
[x1,y1]=ginput;  %Beenden durch Enter-Taste
%Es werden Punkte durch Anklicken auf [x1 y1]
%\"ubertragen: (x1(j),y1(j)), j=1,...
%Das  erste x1 (=x1(1)) wird weggelassen und durch das 
%zweite x1 (=x1(2)) ersetzt, etc. Die y1-Werte werden ignoriert.
Bilddatei_erzeugen=['print -dps Fall',num2str(Fall),'Grad',num2str(n-1),'.ps'];
Punktzahl=length(x1);  %mu\ss\ gerade sein
for j=1:Punktzahl/2
  if x1(j)<a | x1(j)>b, %geklickte Werte au\ss erhalb [a,b], dann aufh\"oren
    Differenz=oben-unten;Differenzrel=2*Differenz/(oben+unten);
    disp(['Differenz (relativ zum Fehler) zwischen oberer und unterer ',...
          'Schranke (nicht der Fehler): ',num2str(Differenzrel)]);
    disp(['Fehler ca. ',num2str(0.5*(oben+unten)),', Rechnung beendet']);
    hold off;
    %disp(['Eine Bilddatei wird mit ',Bilddatei_erzeugen,' erzeugt!']);
    %eval(Bilddatei_erzeugen) %Es wird eine Postskript-Datei vom Typ
                              %"FalljGradn.ps" erzeugt
    return;
  end; %if
  [minwert,kn]=min(abs(Knoten-x1(2*j-1))); %nicht so emfindlich
  %kn=find(abs(Knoten-x1(2*j-1))<=0.05*x1(2*j-1)); %empfindlich!
  Knoten(kn)=x1(2*j);
  disp(['Ersetzt wird Knoten Nr. ',num2str(kn),' durch ',num2str(Knoten(kn))]);
  y=remeztestf(Fall,Knoten(kn));
  yy=polyval(p,Knoten(kn));
  Feh=y-yy;
  plot([Knoten(kn),Knoten(kn)],[0,Feh],'--');
end; %for j
hold off 
Differenz=oben-unten;Differenzrel=2*Differenz/(oben+unten);
    disp(['Differenz (relativ zum Fehler) zwischen oberer und unterer ',...
          'Schranke (nicht der Fehler): ',num2str(Differenzrel)]);
    disp(['Fehler ca. ',num2str(0.5*(oben+unten))]);
return;
