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

1.2. Метод Симпсона для интегрирования функции от двух переменных F(x,y) по прямоугольной двумерной области и его реализация на языке Pascal.

Функция double_simpson (листинг 1.3.) является, по существу, прямым обобщением одномерного метода Симпсона на случай вычисления двойного интеграла по прямоугольной области.

Листинг 1.3. Функции double_simpson и simple_simpson модуля integral
Function double_simpson(F:real_fun2; x0,x1,y0,y1:real; x_div,y_div:word):real;
var
dx,dy,x,sum:real;
i:word;
function simple_simpson(x:real):real;
var
y,sum:real;
j,v:word;
begin
sum:=F(x,y0)+F(x,y1);

y:=y0;
for j:=1 to 2*y_div-1 do

begin
y:=y+dy;
if Odd(j) then
sum:=sum+4.0*F(x,y)
else
sum:=sum+2.0*F(x,y);
end;
simple_simpson:=sum;
end;{simple_simpson}
begin{doudle_simpson}
dx:=(x1-x0)/(2.0*x_div);
dy:=(y1-y0)/(2.0*y_div);
x:=x0;
sum:=simple_simpson(x0)+simple_simpon(x1);
for i:=1 to 2*x_div-1 do
begin
x:=x+dx;
if Odd(i) then
sum:=sum+4.0*simple_simpson(x)
else
sum:=sum+2.0*simple_simpson(x);
end;
double_simpson:=dx+dy*sum/9.0;
end;{double_simpson}

Недостатком рассмотренных функций интегрирования является то, что они не дают возможности явно задать точность вычисления интеграла. Точность связана с количеством точек разбиения, но ее значение в этих функциях не определяется с адаптивным выбором шага разбиения. Такой функцией является adaptive_simpson. Параметры eps и eta задают соответственно абсолютную и относительную погрешности. Их роль поясняется следующим неравенством:

.
Функция adaptive_simpson (листинг1.4) использует рекурсивную процедуру simpson3point, которая вычисляет значение интеграла по интервалу [x₀, x₀+δx], где x₀ – не обязательно исходная левая граничная точка.
Если трехточечный метод Симпсона не дает достаточную точность на данном интервале, этот интервал делится на три равные части, и метод вновь применяется к каждой из полученных частей. В результате получим 7 точек разбиения, но вычислять функцию F(x) придется только в четырех из них, поскольку значения в других трех точках уже известны.
При адаптивном разбиении имеется одна тонкость. При переходе к подынтервалам, составляющим одну треть от исходного, чтобы получить новые абсолютную и относительную погрешности, надо поделить eps и eta на .