%PROGRAM Romberg;
%  USES Crt;
%  CONST  M = 10;
%  TYPE Matrix = ARRAY[0..M,0..M] OF Real;
%  VAR T: Matrix;
%  Var N, NN, i, ii, k, p : Integer;
%                   Count : Integer;
%                 Pi, Pi2 : Real;
%          h, a, b, x, eps: Real;
%                      tt : text;
%
%   FUNCTION Abbruch : Boolean;
%   VAR ch: Char;
%   BEGIN
%     Write('Abbruch ja oder nein ?? ');
%     ch := ReadKey;
%     Abbruch := (ch = 'J') OR (ch = 'j')
%   END; { Abbruch }
%
%   FUNCTION f(x: Real): Real;
%     BEGIN
%       f := 1;
%     END;  (* function *)
% 
%   FUNCTION g(x: Real): Real;
%     BEGIN
%       (*g := 2 * Exp(x*x);
%       g := Sin(x) / Cos(x);  *)
%       (* Berechnung des Umfangs einer Ellipse *)
%       (*g := Sqrt(1-3/4*Sin(Pi2*x)*Sin(Pi2*x));*)
%       g:=1/(1+sin(x)*sin(x));
%       Count:= Count+1;
%     END;  (* function *)

%  BEGIN  (* main *)
% 
%       Pi := 4* Arctan(1.0);
%       Pi2 := Pi/2.0;
%       Count := 0;
% 
%       Assign(tt,'m:\user\pascal\dat\romberg.out');
%       Rewrite(tt);                       (* neue Datei er”ffnen *)
% 
%     REPEAT
%       Writeln;
% 
%       FOR k:=0 TO M DO
%         FOR i:=0 TO M DO T[k,i] := 0.0;
% 
%       N := 0;      (* Anzahl Halbierungen *)
%       k := 0;
%       NN := 1;
%       Write('a = ');
%       Readln(a);
%       Write('b = ');
%       Readln(b);
%       Write('Toleranz eps =  ');
%       Readln(eps);
%       h := (b-a);  (* Intervall-Laenge *)
%       T[0,0] := h*(g(a) + g(b))/2.0;
%       Writeln(tt, 'Kontroll-Output f?r Romberg:');
%       Writeln(tt, 'a = ', a, ' b = ', b, ' b-a = ', h);
%       Writeln(tt, 'eps = ', eps);
%       Writeln(tt,Pi);
% 
%       REPEAT
%         N := N + 1;
%         k := k + 1;
%         T[k,0] := T[k-1,0]/2.0;
%         h := h/2.0;
%         NN := NN*2;
%         FOR i := 1 TO (NN DIV 2) DO              (* NN/2 = Anz. neue Mittelpkte *)
%         BEGIN
%           ii := 2*i - 1;
%           T[k,0] := T[k,0] + h*g(a+ii*h)
%         END;
%         (* Romberg *)
%         p := 1;                 (* Hilfsgroesse fuer 4k *)
%         FOR i := 1 TO k DO
%         BEGIN
%           p := 4*p;
%           T[k,i] := (p*T[k,i-1] - T[k-1,i-1])/(p - 1)
%         END;
% 
%       UNTIL Abs(T[k,k-1] - T[k-1,k-1]) < eps;
%       FOR k := 0 TO M DO
%       BEGIN
%         Writeln(tt);
%         FOR i := 0 TO k DO Write(tt, T[k,i]:13:12,' ')
%         (* Faktor 4*Pi f?r den Umfang der Ellipse *)
%         (*FOR i := 0 TO 3 DO Write(tt, 4*Pi*T[k,i]:13:12,' ') *)
%       END;
%       Writeln(tt);
%       Writeln(tt, 'Anzahl Funktionsauswertungen: ', Count);
% 
%     UNTIL Abbruch;
% 
%     Close(tt)                         (* Datei permanent speichern
%                                        im File romberg.out
%                                        im directory m:\user\pascal\dat *)
%
%  END. (* main *)


% Beispiel ... Serie 9 SoSe 2004
clear all
hold off

bsp_nummer = input('bsp_nummer (1 2 3 oder 4) = ')

switch bsp_nummer 
    case 1
     f = inline('1/x')                   % a = 1 und b = 2
     a = 1
     b = 2
     int = log(2)
    case 2
     f = inline('sin(x)/(1 + x^2)')      % a = 0 und b = pi/2
     a = 0
     b = pi/2
     int = .52697855761398923223;
    case 3
     f = inline('exp(x)*cos(x)')         % a = 0 und b = pi/2
     a = 0
     b = pi/2
     int = 1/2*(exp(pi/2) - 1)
    case 4
     f = inline('sqrt(x)')         % a = 0 und b = 1
     a = 0
     b = 1
     int = 2/3

end

% a = input('untere Intervallgrenze a = ');
% b = input('obere  Intervallgrenze b = ');
eps = input('Toleranz eps = ');
Nmax = input('maximale Anzahl Halbierungen Nmax = ');

k = 2;

h = b - a;              % h0
T0 = h/2*(f(a) + f(b));
M0 = h*f(a + h/2);
T1 = 1/2*(T0 + M0);
h = h/2;                % h1
n = 2;
R(1,1) = T0;
R(2,1) = T1;
R(2,2) = 1/3*(4*R(2,1) - R(1,1));
quot = [];
while abs(R(k,k) - R(k-1,k-1)) > eps & k < Nmax
    k = k + 1;
    h = h/2;            % h2
    n = 2*n;
    mitten = 0.0;
    for i = 1:2:n
         mitten = mitten + f(a+i*h);
    end
    mitten = 2*h*mitten;
    R(k,1) = 1/2*(R(k-1,1) + mitten);
    
    % Romberg Schema
    p = 1;
    for j = 2:k
        p = 4*p;
        R(k,j) = (p*R(k,j-1) - R(k-1,j-1))/(p-1);
        
    end
    
   
end


 % exakter Wert (Maple)
 % Berechnung der Quotienten q
 % vgl. Roy Black
 

    mm = length(diag(R))
    
    for k = 2:mm
        q(k) = (R(k,k) - int)/(R(k-1,k-1) - int);
    end
    figure(1)
    %hold off
    semilogy(2:mm, abs(q(2:mm)), 'r.-')
    title('Verfahren von Romberg mit fortgesetzter Halbierung', 'interpreter', 'latex')
    xlabel('Anzahl Halbierungen', 'interpreter', 'latex')
    ylabel('Konvergenz der Quotienten $|q_{k}|$', 'interpreter', 'latex')
    grid on
