%
% -----------------------------------------------------
%
%   Approximation                   Wintersemester 2013/14
%   (H.J. Oberle)
%   T-Bestapproximation mittes des REMEZ-Algorithmus
%
%   Aufgabe 41 a)
%   f(t) = 1/(1+t), 0 <= t <= 1,
%   p(t) = a0 + a1*t + a2*t^2
%
% -----------------------------------------------------
%
function aufg41a
%
clear all
close all
format long e
format compact
%   
ta  = 0;
tb  = 1;
x   = [0 0.25 0.75 1];
%   
n1  = length(x);
n   = n1 - 1;
%
app = zeros(n,1); 
itmax = 20;
tol   =  1.e-8;  
%
[x,app,mue] = remez(x,ta,tb,itmax,tol);
%
t = [ta:0.001:tb];
y = fcn(t);
p = approx(t,app);
e = y-p;
%
plot(t,e,'b','Linewidth',1.2)
set(gca,'FontSize',12)
xlabel('Bestapproximation; Aufg.41 a)','FontSize',14)
ylabel('Fehlerfunktion e(t)','FontSize',14)
axis([0 tb -8e-3 8e-3])
hold
for i=1:n1
    z(i) = (-1)^(i-1)*mue;
end
plot(x,z,'b.','Markersize',20) 
plot([ta tb],[0 0],'--','Linewidth',1.2)
%
print -depsc aufg41a
%
%
% -------------------------------------------------------------------
%
function [y] = fcn(t)
%
%   Die zu approximierende Funktion
%   t :   Abszisse
%   y :   Funktionswert in x 
%
    y    =  (1+0*t)./(1+t);
%
return;
%
% -------------------------------------------------------------------
%
function [H0] = basis(x)
%
%   Die  approximierende Funktion
%   x :   Referenz
%   H0:   Auswertung der Basis auf der Referenz
%
    H0   =  [(1+0*x); x; (x.*x)]';                     %'
%
return;
% -------------------------------------------------------------------
%
function [p] = approx(x,app)
%
%   Die  approximierende Funktion
%   x :   Abszisse
%   p :   Funktionswert in x 
%   app : Koeffizienten der Basisdarstellung
%
    p    =  ((app(3)*x+app(2)).*x)+app(1);
%
return;
%
% -------------------------------------------------------------------
%
function [x,app,mue] = remez(x,ta,tb,itmax,tol)
%
n1   =  length(x);
n    =  n1 - 1;
iter =  1;
konv =  0;
v    =  zeros(n1,1); 
for i=1:n1
    v(i) = -(-1)^i;
end
fprintf('\n   Toleranz:  %3.5e\n',tol)
%
while and(konv==0,iter < itmax)
%
   fprintf('\n\n   -----   Iteration:  %3.0f  -----------------\n',iter)
   fprintf(' Referenz: ')
   fprintf('\n  %+3.5e %+3.5e %+3.5e %+3.5e %+3.5e',x)
%  
   H0 =  basis(x);   
   y  =  fcn(x); 
%
%  --- Loesung des linearen Gleichungssystems
%
   AM   =  [H0 v];
   d    =  AM\y';                                         
%
   app  =  d(1:n);
   mue  =  d(n1);   
   fprintf('\n Referenzloesung: ')
   a = d';                                                      
   fprintf('\n  %+3.5e %+3.5e %+3.5e %+3.5e %+3.5e',a)
   ip   =  1;
   if mue < 0
      ip = 0;
   end
%
%  --- Bestimmung der Maxima/Minima in jedem Teilintervall
%
   t0 = ta;
   [y0] = fcn(t0);
   [p0] = approx(t0,app);
   e0 = y0-p0;
%
   tmax  = t0;
   emax  = e0;
   kint   = 0;
   if t0 >= x(1)
      kint = 1;
   end
   kmax  =  kint;
%
   intmax = n;
   if x(n1) < tb
      intmax = n1;
   end
%
%  --- Die While-Schleife durchlaeuft die Intervalle
%
   while kint <= intmax,  
%
      if kint < n1
         t1 = x(kint+1);
      else
         t1  = tb;
      end
%
%  --- Absuchen auf feinem Gitter
%
      h = (t1 - t0)*0.001;
      t = [t0:h:t1];
      y = fcn(t);
      p = approx(t,app);
      e = y-p;
      [ee,k] = max(abs(e));
      tk  =  t(k);
      ek  =  e(k);
%
%  --- Festhalten des maximalen Fehlerwertes:
%
      if abs(ek) > abs(emax)
         tmax  = tk;
         emax  = ek;
         kmax  = kint;
      end
%
%  --- Ende der while - Schleife:
      t0  =  t1;
      kint = kint + 1;
%   
   end % while kint
%
   fprintf('\n  Maximaler Fehler:  tmax, emax, kmax: \n')
   fprintf('   %+2.5e   %+2.5e  %6.0f \n', tmax, emax, kmax)  
%
%  ---  Konvergenzkriterium
%
   if abs(emax) - abs(mue) < tol*abs(mue)
      konv  =  1;
      fprintf('\n\n  Konvergenzkriterium ist erfuellt! \n')
   end 
%
%  --- Austausch eines Referenzpunktes:
%
   if kmax == 0           
     if or(and(ip==1,emax>0),and(ip==0,emax<0))
         x(1) = tmax;
     else
       for i=1:n
          x(n1+1-i) = x(n1-i);
       end
       x(1) = tmax;
     end 
   end
%
   if and(kmax>0,tmax <= x(n1))        
     if and((-1)^(kmax+ip)>0,emax<0)
        x(kmax+1) = tmax;
      end
     if and((-1)^(kmax+ip)>0,emax>0)
        x(kmax) = tmax;
     end
     if and((-1)^(kmax+ip)<0,emax<0)
        x(kmax) = tmax;
     end
     if and((-1)^(kmax+ip)<0,emax>0)
        x(kmax+1) = tmax;
     end
   end
%
   if tmax > x(n1)   
     if or(and((-1)^(ip+n1)<0,emax<0),and((-1)^(ip+n1)>0,emax>0))
         x(n1) = tmax;
     else
       for i=1:n
          x(i) = x(i+1);
       end
       x(n1) = tmax;
     end
   end
%
%
   iter  =  iter + 1;
%
end % while - Iteration
%
return;
% -------------------------------------------------------------------
%





