Romberg-Integration

program Romberg(input, output);

    {Romberg-Integration der Funktion f(x) = exp(sin^2(x))}
    
    const
        MaxArray = 50;
        
    var
        ints:       array [1..MaxArray, 1..MaxArray] of double;
        a, b:       double;
        eps, epsa:  double;
        maxit:		integer;
        n:          longint;
        i, j, k:	longint;
        
    function f(x: double): double;
    begin {f}
        f := exp(sqr(sin(x)))
    end; {f}
            
    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
        	writeln('Romberg-Integration der Funktion f(x) = exp(sin^2(x))');
        	writeln;
        	write('Wie lautet die erste Intervallgrenze? ');
        	readln(a);
        	write('Wie lautet die zweite Intervallgrenze? ');
        	readln(b);
        	writeln;
        	write('Wie gross ist das Abbruchkriterium? ');
        	readln(eps);
        	write('Wie lautet die maximale Iterationstiefe? ');
        	readln(maxit);
        	if maxit > MaxArray then begin
        	    writeln('Die Iterationstiefe ist zu gross (maximal ', MaxArray,
        	        ').');
        	    halt
        	end; {if}
        	writeln
         end; {init}

    procedure result;
        begin {result}
            writeln('Das Integral hat den Wert ', ints[1, i + 1]:7:3)
        end; {result}

    function Trapezregel(n: longint; a, b: double): double;
        var
            sum, x: double;
            step:   double;
            i:      integer;
        begin {Trapezregel}
            sum := 0;
            step := (b - a) / n;
            for i := 1 to n - 1 do begin
                x := a + i * step;
                sum := sum + f(x)
            end; {for}
            Trapezregel := (b - a) * (f(a) + 2 * sum + f(b)) / (2 * n)
        end; {Trapezregel}

    function potenz(a, n: integer): longint;
        begin {potenz}
            potenz := round(exp(n * ln(2)))
        end; {potenz}
        
    begin {Romberg}
        init;
        ints[1, 1] := Trapezregel(1, a, b);
        epsa := 1.1 * eps;
        i := 0;
        while (epsa > eps) and (i < maxit) do begin
            i := i + 1;
            n := potenz(2, i);
            ints[i + 1, 1] := Trapezregel(n, a, b);
            for k := 2 to i + 1 do begin
                j := 2 + i - k;
                ints[j, k] := potenz(4, k - 1) * ints[j + 1, k - 1] - ints[j, k - 1];
                ints[j, k] := ints[j, k] / (potenz(4, k - 1) - 1)
            end; {for}
            epsa := abs((ints[1, i + 1] - ints[1, i]) / ints[1, i + 1]) * 100
        end; {while}
        result     
    end. {Romberg}