Приближенное вычисление интеграла

Листинг 1.5. Функция romberg модуля integral
Function romberg (F:real_fun;x0,x1,eps,eta:real; min,max:word):real;
const
abs_max=30;
var
p,dx,error,F_of_x0, F_of_x1, F_of_xk,
roundoff_error,integral_abs,tolerance,
previous_estimate,current_estimate,
mid_sum, temp_sum, mid_sum_abs:real;
table:array[0..abs_max] of real;
j,n:word;
k,r:longint;
done:Boolean;
denom:array[1..abs_max] of real;
begin
p:=1.0;
for k:=1 to abs_max do
begin
p:=4.0*p;
denom[k]:=1.0/(p-1.0);
end;
dx:=x1-x0;
F_of_x0:=F(x0);
F_of_x1:=F(x1);
current_estimate:=0.0;
previos_estimate:=0.0;
done:=False;
table[0]:=0.5*dx*(F_of_x0+F_of_x1);
integral_abs:=0.5*Abs(dx)*(Abs(F_of_x0)+Abs(F_of_x1));
n:=1;
r:=1;
repeat
dx:=0.5*dx;
mid_sum:=0.0;
mid_sum_abs:=0.0;
roundoff_error:=0.0;
for k:=1 to r do
begin
F_of_xk:=F(x0+(2*k-1)*dx);
mid_sum_abs:=mid_sum_abs+Abs(F_of_xk);
F_of_xk:=F_of_xk+roundof_eroor;
temp_sum:=mid_sum+F_of_xk;
roundof_error:=(mid_sum-temp_sum)+F_of_xk;
mid_sum:=temp_sum;
if KeyPressed then
Halt;
end;
table[n]:=0.5*table[n-1]+dx*mid_sum;
integral_abs:=0.5*integral_abs+Abs(dx)*mid_sum_abs;
for j:=n-1 downto 0 do
table[j]:=table[j+1]+denom[n-j]*(table[j+1]-table[j]);
if n>=min then
begin
tolerance:=eta*integral_abs+eps;
error:=Abs(table[0]-current_estimate)+Abs(current_estimate-previos_estimate);
done:=(errorend;
Inc(n);
done:=done or(n>max);
previous_estimate:=current_estimate;
current_estimate:=table[0];
r:=r+r;
until done;
romberg:=current_estimate;
end;


1.4. Метод Гаусса и его реализация на языке Pascal.

Теперь перейдем к гауссовским квадратурам – семейству правил интегрирования, основанных на неравномерном разбиении основного интервала интегрирования. Вообще, метод гаусса с n точками точен для полиномов степени 2n – 1. В функции gauss3 (листинг1.6.) основной трехточечный алгоритм Гаусса применяется к каждой из n равных частей интервала. Для интервала [-1,1] узлами квадратурной формулы являются нули полинома Лежандра третьей степени P₃ = (5x³ – 3x)/2, а коэффициенты выбираются специальным образом.