Вычисление интеграла

const
  A = 1;
  B = 81;
  Res = 1.595;
 
function F(X: Real): Real;
begin
  F := 1/(4*Sqrt(X)*(3+Exp(1/4)*Ln(X)));
end;
var
  I, N: Integer;
  S, X, H, Eps, D, W: Real;
begin
  WriteLn('Введите погрешность:');
  ReadLn(W);
  D := 1;
  while ((D > W) or (N < 1000)) do
    begin
      X := A;
      N := N * N;
      H := (B-A)/N;
      S := 0;
      D := 0;
      Eps := 0;
      for I:=1 to N do
        begin
          S := S + F(X) * H;
          X := X + H;
        end;
    D := Abs(Res - S);
    Eps := Abs((D / Res) * 100);
  end;
  WriteLn('По формуле прямоугольников с недостатком = ', S:8:4);
  WriteLn('Абсолютная погрешность достигнута при количестве узлов = ', N);
  WriteLn('Абсолютная погрешность = ', D:10:8);
  WriteLn('Относительная погрешность = ', Eps:8:4); 
  N := 0;
  D := 1;
  while D > W do
    begin
      X := A;
      N := N + 1;
      H := (B-A)/N;
      S := 0;
      D := 0;
      Eps := 0;
      for I:=1 to N do
        begin
          S := S + F(X) * H;
          X := X + H;
        end;
    D := Abs(Res - S);
    Eps := Abs((D / Res) * 100);
  end;
  WriteLn('По формуле прямоугольников с избытком = ', S:8:4);
  WriteLn('Абсолютная погрешность достигнута при количестве узлов = ', N);
  WriteLn('Абсолютная погрешность = ', D:10:8);
  WriteLn('Относительная погрешность = ', Eps:8:4);
  N := 0;
  D := 1;
  while D > W do
    begin
      X := A;
      N := N + 1;
      H := (B-A)/N;
      S := 0;
      D := 0;
      Eps := 0;
      for I:=1 to N do
        begin
          S := S + F(X) * H;
          X := X + H;
        end;
    D := Abs(Res - S);
    Eps := Abs((D / Res) * 100);
  end;
  WriteLn('По формуле трапеций = ', S:8:4);
  WriteLn('Абсолютная погрешность достигнута при количестве узлов = ', N);
  WriteLn('Абсолютная погрешность = ', D:10:8);
  WriteLn('Относительная погрешность = ', Eps:8:4);
  N := 0;
  D := 1;
  while D > W do
    begin
      X := A;
      N := N + 1;
      H := (B-A)/N;
      S := F(A) + F(B);
      D := 0;
      Eps := 0;
      for I:=1 to N-1 do
        begin
          X := X + H;
          if (I mod 2) = 0 then
            S := S + (F(X) * 2)
          else
            S := S + (F(X) * 4)
        end;
      S := S * H/3;  
      D := Abs(Res - S);
      Eps := abs((D / Res) * 100);
    end;
  WriteLn('По формуле Симпсона = ', S:8:4);
  WriteLn('Абсолютная погрешность достигнута при количестве узлов = ', N);
  WriteLn('Абсолютная погрешность = ', D:10:8);
  WriteLn('Относительная погрешность = ', Eps:8:4);
end.