Вычислить определенный интеграл по формуле Ньютона-Котеса - Pascal

Формулировка задачи:

Здравствуйте! Надо вычислить определенный интеграл по формуле Ньютона-Котеса открытого типа с заданной точностью (eps =10^(-4)) c использованием АДАПТИВНОГО АЛГОРИТМА! Перерыла весь интернет и Ваш форум. Чисто логически поняла, как это сделать, но вот с написанием кода проблемы(( Подскажите, пожалуйста, как нужно изменить прогу? Вот что получилось (без адаптивного алгоритма)

program N-K;
const n = 8; //Кол-во разбиений
a = 0; //начало
d = 1; //конец
h = 1/9; //шаг
e = power(10,-4); //точность
 
c : array of real = (1787/4480, -2803/4480, 4967/4480, -1711/4480, -1711/4480, 4967/4480, -2803/4480, 1787/4480); //Коэф 
var
i : integer;
R, x : real;
 
//Вычисление значения функции
function f(x : real) : real;
begin
f := ln(1+x)/(sqr(x)+1);
end;
 
begin
x :=0;
R :=0;
for i :=0 to n-1 do
begin
x := x+h;
R := R + c[i]*f(x);
end;
 
writeln(R); 
 
readln;
end.

Решение задачи: «Вычислить определенный интеграл по формуле Ньютона-Котеса»

Листинг программы

{
 Программа для вычисления интеграла
 численным методом Ньютона-Котеса открытого типа 4-го порядка
 с заданной погрешностью вычислений
}
program NewtonCotes;
 
  {подинтегральная функция}
  function F(x: real): real;
  begin
    F := sin(x);
  end;
 
  {вычисление интеграла на отрезке [a, b] методом Ньютона-Котеса открытого типа 4-го порядка}
  function NewtonCotesOpen4(a, b: real): real;
  const
    n = 4; {порядок метода}
    c: array [1..n - 1] of real = (2, -1, 2); {весовые коэффициенты}
    K = 3.0; {делитель}
  var
    Res: real;
    h: real; {величина шага}
    x: real;
    i: integer;
  begin
    h := (b - a) / n;
    Res := 0;
    for i := 1 to n - 1 do
    begin
      x := a + i * h;
      Res := Res + c[i] * F(x);
    end;
    NewtonCotesOpen4 := Res * (b - a) / K;
  end;
 
  {вычисление интеграла на отрезке [a, b] как суммы интегралов на n подинтервалах}
  function SumSubIntervals(a, b: real; n: integer): real;
  var
    h: real;
    i: integer;
  begin
    SumSubIntervals := 0;
    h := (b - a) / n;
    for i := 1 to n do
      SumSubIntervals := SumSubIntervals + NewtonCotesOpen4(a + pred(i) * h, a + i * h);
  end;
 
var
  a, b,        {границы интегрирования}
  Eps: real;   {погрешность численного интегрирования}
  n: integer;  {количество подинтервалов интегрирования}
  Icur,        {значение интеграла на текущем шаге итерации}
  Ipre: real;  {значение интеграла на предыдущем шаге итерации}
 
begin
  a := 0;
  b := Pi;
  Eps := 0.00001;
 
  n := 1;
  Icur := SumSubIntervals(a, b, n);
  repeat
    Ipre := Icur;
    n := 2 * n;                {увеличиваем количество подинтервалов}
    Icur := SumSubIntervals(a, b, n);
    WriteLn('I=', Icur: 10: 8, ' при n=', n, ', e=', Abs(Icur - Ipre): 10: 8);
  until Abs(Icur - Ipre) < Eps;
  writeln('------------------------------------------------------');
  WriteLn('I=', Icur: 10: 8, ' при n=', n, ', e=', Abs(Icur - Ipre): 10: 8);
end.

Объяснение кода листинга программы

{Программа для вычисления интеграла численным методом Ньютона-Котеса открытого типа 4-го порядка} program NewtonCotes; {подинтегральная функция} function F(x: real): real; begin F := sin(x); end; {вычисление интеграла на отрезке [a, b] методом Ньютона-Котеса открытого типа 4-го порядка} function NewtonCotesOpen4(a, b: real): real; const n = 4; {порядок метода} c: array [1..n - 1] of real = (2, -1, 2); {весовые коэффициенты} K = 3.0; {делитель} var Res: real; h: real; {величина шага} x: real; i: integer; begin h := (b - a) / n; Res := 0; for i := 1 to n - 1 do begin x := a + i h; Res := Res + c[i] F(x); end; NewtonCotesOpen4 := Res (b - a) / K; end; {вычисление интеграла на отрезке [a, b] как суммы интегралов на n подинтервалах} function SumSubIntervals(a, b: real; n: integer): real; var h: real; i: integer; begin SumSubIntervals := 0; h := (b - a) / n; for i := 1 to n do begin SumSubIntervals := SumSubIntervals + NewtonCotesOpen4(a + pred(i) h, a + i h); end; end; var a, b, Eps: real; n: integer; Icur, Ipre: real; begin a := 0; b := Pi; Eps := 0.00001; n := 1; Icur := SumSubIntervals(a, b, n); repeat Ipre := Icur; n := 2 n; Icur := SumSubIntervals(a, b, n); WriteLn('I=', Icur: 10: 8, ' при n=', n, ', e=', Abs(Icur - Ipre): 10: 8); until Abs(Icur - Ipre) < Eps; writeln('------------------------------------------------------'); WriteLn('I=', Icur: 10: 8, ' при n=', n, ', e=', Abs(Icur - Ipre): 10: 8); end.