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

Листинг 1.7. Процедура compute_gauss_coeffs модуля integral
procedure compute_gauss_coeffs(deg:word);
const
eps=6.0e-20;
var
i,index:word;
P0k, P0k_1,D0k, P1k,P1k_1, D1k,
x0,x1,y,z,dx,x,u:real;
procedure legendre_poly(n:word; x:real; var Pk,Pk_1,Dk:real);
var
Pk_2,Dk_1,Dk_2:real;
i,j,k:word;
begin
If n=0 then
Begin
Pk:=1.0;
Dk:=0.0;
end
else
begin
Pk_1:=1.0;
Pk:=x;
Dk_1:=0.0;
Dk:=1.0;
i:==3;
j:=1;
for k:=2 to n do
begin
Pk_2:=Pk_1;
Pk_1:=Pk;
Dk_2:=Dk_1;
Dk_1:=Dk;
Pk:=(i*x*Pk_1-j*Pk_2)/k;
Dk:=(i*(Pk_1+x*Dk_1)-j*Dk_2)/k;
Inc(I,2);
Inc(j);
end;
end
end;{legendre_poly}
begin{computr_gauss_coeffs}
index:=(deg+1) div 2;
dx:=1.0/(10.0*deg);
x0:=0.0;
x1:=x0+dx;
if Odd(deg) then
begin
zero[index]:=0.0;
legendre_poly(deg,x0,P0k, P0k_1,D0k);
weight[index]:=2.0/(P0k_1*D0k*deg);
end;
for i:=0 to 10*deg-1 do
begin
x0:=x1;
x1:=x1+dx;
legendre_poly(deg,x0,P0k_1,D0k);
legendre_poly(deg,x1,P1k,P1k_1,D1k);
if P0k*P1k<=0.0 then
begin
x:=x0-P0k*dx/(P1k-P0k);
legendre_poly(deg,x,P0k,P0k_1,D0k);
u:=P0k/D0k;
y:=x-u;
while Abs(x-y)>=eps do
begin
if keyPressed then begin
writeln(‘>=eps loop:’,x:10:10,’ ‘, y:10:10,’ ‘,Abs(x-y):10);
readln;
end;
x:=y;
legendre_poly(deg,x,P0k,P0k_1,D0k);
u:=P0k/D0k;
y:=x-u;
end;
inc(index);
legendre_poly(deg,y,P0k,P0k_1,D0k);
zero[index]:=y;
weight[index]:=2.0/(P0k_1*D0k*deg);
if index=deg then
Break;
end;
end;
For i:=1 to deg div 2 do
begin
Zero[i]:=zero[deg-i+1];
Weight[i]:=weight[deg-i+1];
end;
end;{compute_gauss_coeffs}

В функции gauss (листинг 1.8.) запрограммирован один гауссовский шаг на заданном интервале. Конечно, серьезная прикладная программа будет делить интервал на меньшие подынтервалы, применять эту процедуру к каждому из них адаптивным способом.