%
% -----------------------------------------------------
%
%   Approximation                   Wintersemester 2013/14
%   (H.J. Oberle)
%   T-Bestapproximation mittes des REMEZ-Algorithmus
%
%   Aufgabe 41 b)
%   f(t) = (sin t)/t, 0 <= t <= pi/6,
%   p(t) = a0 + a1*t^2 + a2*t^4 + a3*t^6 
%
% -----------------------------------------------------
%
function aufg41b
%
   clear all
   close all
   format long e
   format compact
%   
   ta  =   0;
   tb  =   pi/6;
   x   =  [0 0.3*tb 0.5*tb 0.8*tb tb];
%   
   n1  =  length(x);
   n   =  n1 - 1;
%
   app = zeros(n,1); 
   itmax = 50;
   tol   =  1.e-6;  
%
   [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 b)','FontSize',14)
   ylabel('Fehlerfunktion e(t)','FontSize',14)
   axis([0 tb -1.5e-10 1.5e-10])
   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 aufg41b
%
%
% -------------------------------------------------------------------
%
function [y] = fcn(t)
%
%   Die zu approximierende Funktion
%   t :   Abszisse
%   y :   Funktionswert in x 
%
    for i=1:length(t)    
    if t(i) < 1.e-5
       y(i)  =  1 - (1/6)*t(i)*t(i);
    else
       y(i)  =  sin(t(i))/t(i);
    end
    end
%
return;
%
% -------------------------------------------------------------------
%
function [H0] = basis(x)
%
%   Die  approximierende Funktion
%   x :   Referenz
%   H0:   Auswertung der Basis auf der Referenz
%
    H0   =  [(1+0*x); x.*x; x.^4; x.^6]';     
%
return;
% -------------------------------------------------------------------
%
function [p] = approx(t,app)
%
%   Die  approximierende Funktion
%   t :   Abszisse
%   p :   Funktionswert in x 
%   app : Koeffizienten der Basisdarstellung
%
    t2   =  t.*t;
    p    =  (((app(4)*t2 + app(3)).*t2+app(2)).*t2) + 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;
% -------------------------------------------------------------------
%






