function A28(Funktion, Verfahren)
% A28(Funktion, Verfahren)
%
% Untersuchung des Konvergenzverhaltens von Quadraturformeln
% für die Funktionen
%              0) f(t) = |t|^3
%              1) f(t) = 1/(1+t^2)
%              2) f(t) = t^10
% und den Verfahren:
%              0) Gauss-Legendre-Quadratur
%              1) Trapez-Regel
%              2) Simpson-Regel
%              3) Romberg-Verfahren

C = ['k','r','g','b'];

for Funktion=0:2
   figure;
      
   if Funktion==0
      Integral = 0.5;
      Func = 'f(t) = |t|^3';
   elseif Funktion==1
      Integral = pi/2;
      Func = 'f(t) = 1/(1+t^2)';
   else
      Integral = 2/11;
      Func = 'f(t) = t^{10}';
   end;
   
   for Verfahren=0:3
      if Verfahren<=1
         nVec = 1:64;
      elseif Verfahren==2
         nVec = 2:2:64;
      else
         nVec = 2.^(0:6)+1;
      end;
      
      err = zeros(size(nVec));
      for Count=1:length(nVec)         
         if Verfahren<3
            n = nVec(Count);
            h = 2/n;
            
            if Verfahren==0
               % b_(k+1)^2 = k^2 / ((2k-1)*(2k+1))
               %           = k^2 / (4k^2-1)
               %           = 1 / (4-1/k^-2)
               
               b     = 1./sqrt(4-1./(1:n-1).^2); % Nebendiagonalen
               T     = diag(b,1)+diag(b,-1);     % Matrix aufsetzen
               [V,D] = eig(T);                   % Eigenwerte und -vektoren berechnen
               x     = diag(D)';                 % Eigenwerte = Stützstellen
               alpha = 2*V(1,:).^2;              % 1. Komponenten der Eigenvektoren = Koeffizienten
               
               if (n==40) & (Funktion==0)
                  % Koeffizienten über Stützstellen berechnen
                  H = gcf;
                  figure;
                  plot(x,alpha,'.k');
                  drawnow;
                  figure(H);
               end;         
            elseif (Verfahren==1) | (Verfahren==2)
               x     = linspace(-1,1,n+1);
               alpha = h*[0.5 ones(1,n-1) 0.5];
               
               if Verfahren==2   % Simpson-Regel mit der Trapez-Regel berechnen
                  % Sh = 1/3 (4 Th - T2h)
                  Th    = alpha;
                  alpha = 4/3*Th;
                  alpha(1:2:end) = alpha(1:2:end) - Th(1:2:end)*2/3;
               end;               
            end;
            
            Quad = sum(alpha.*f(x,Funktion));
         else
            % Romberg-Verfahren
            
            % Erst die Trapezsummen T_h0,...,T_hm ausrechnen...
            m    = log2(nVec(Count)-1);
            T    = zeros(m+1,1);
                        
            for j=0:m
               n      = 2^j+1;
               h      = 2/(n-1);
               x      = linspace(-1,1,n);
               alpha  = h*[0.5 ones(1,n-2) 0.5];
               T(j+1) = sum(alpha.*f(x,Funktion));
            end;
            
            % ... dann die Grenzwertextrapolation
            for j=1:m
               T(1:m+1-j) = (4^j*T(2:m+2-j) - T(1:m+1-j)) / (4^j-1);
            end;            
               
            Quad = T(1);
         end;
         
         err(Count) = abs(Quad-Integral);
      end;
      
      semilogy(nVec,err+1e-20,strcat('.-',C(Verfahren+1)));
      hold on;
      drawnow;
   end;
   
   legend('Gauss-Legendre','Trapez-Regel','Simpson-Regel','Romberg-Verfahren');
   axis([nVec(1) nVec(end) 1e-16 1]);      
   title(strcat(Func));
   xlabel n;
   ylabel Interpolationsfehler;
end;
   
   
   
   
   

function ff = f(t,Funktion)
if Funktion==0
   ff = abs(t).^3;
elseif Funktion==1
   ff = 1./(1+t.^2);
else
   ff = t.^10;
end;


