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

Листинг 1.1. Интерфейсная секция модуля integral.
Unit integral;
Interface
Const
Max_dim =10;
Max_deg=96;
Type
Real_fun=function(x:real):real;
Real_fun2=function(x,y:real):real;
Real_vec=array[1..max_dim+1] of real;
Index=array[1..max_dim+1] of word;
Vec_fun=function(j:word; x:real_vec):real;
Var
no_evaluations, highest_level:word;
function simpson(F:real_fun; x0,x1:real; div_no:word):real;
function double_simpson(F:real_fun2; x0,x1,y0,y1:real; x_div,y_div:word):real;
function adaptive_simpson(F:real_fun;x0,x1,eps,eta:real):real;
function romberg(f:real_fun; x0,x1,eps,eta:real; min,max:word):real;
function gauss3(F:real_fun;x0,x1:real; n:word):real;
procedure compute_gauss_coeffs(deg:word);
function gauss(Freal_fun:x0,x1:real; deg:word):real;


1.1. Метод Симпсона для интегрирования функций F(x) по заданному промежутку и его реализация на языке Pascal.

Перейдем к секции реализации. Она начинается описанием функции simpson. Стоит сказать несколько слов о выборе узлов и коэффициентов квадратурной формулы Симпсона. Идея трехточечного метода Симпсона заключается в следующем.

Пусть x_m – это средняя точка интервала [x₀, x₁] и пусть Q(x) – единственный полином второй степени, который интерполирует (приближает) подынтегральную функцию F(x) по точкам x₀, x_m и x₁. Искомый интеграл аппроксимируется интегралом от функции Q(x):
I≈ .
Это оценка точна, если F(x) является полиномом степени 3.
В функции simpson интервал интегрирования делится на div_no равных частей, а трехточечная формула Симпсона применяется к каждому такому интервалу. Параметрами функции simpson (листинг1.2) являются, по порядку, подынтегральная функция, нижняя и верхняя границы интервала интегрирования и количество подынтервалов.


Листинг 1.2. Функция simpson модуля integral
Implementation
Uses crt;
Var
zero, weight:array [1..max_deg] of real;
Function simpson(F:real_fun; x0,x1:real; div_no:word):real;
Var
x, dx, sum:real;
j:word;
begin
dx:=(x1-x0)/(2.0*div_no);
sum:=F(x0)+F(x1);
x:=x0;
for j:=1 to 2*div_no-1 do
begin
x:=x+dx;
if Odd (j) then
sum:=sum+4.0*F(x)
else
sum:=sum+2.0*F(x);
end;
simpson:=dx*sum/3.0;
end;