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

Листинг 1.4. Функция adaptive_simpson модуля integral
Function adaptive_simpson(F:real_fun:x0,x1,eps,eta:real):real;
const
max_level=35;
var
k,nest_level:word;
integral_abs:real;
function simpson3poin(x0,delta_x, estimate, integral_abs,
eps,eta,left,middle,right:real):real;
var
dx3,sum,eps3,eta3,factor,left_integ,
middle_integ, right_integ,F1,F2,F4,F5:real;
begin
Inc(nest_level);
dx3:=delta_x/3.0;
F1:=F(x0+0.5*dx3);
F2:=F(x0+dx3);
F:=F(x0+2.0*dx3);
F5:=F(x0+2.5*dx3);
Inc(no_evaluations,4);
factor:=dx3/6.0;
left_integ:=factor*(left+4.0*F1+F2);
middle_integ:=factor*(F2+4.0*middle+F4);
right_integ:=factor*(F4+4.0*F5+right);
sum:=left_integ+middle_integ+right_integ;
integral_abs:=integral_abs- Abs(eastimate)+Abs(left_integ)+Abs(middle_integ)+Abs(right_integ);
if (nest_level>1) and ((nest_level=max_level) or
(Abs(sum- estimate)<=eps+eta*integral_abs)) then simpson3point:=sum
else
Begin
If nest_level>highest_level then
Inc(highest_level);
Eps3:=0.577*eps;
Eta3:=0.577*eta;
Left_integ:=simpson3point(x0,dx3,left_integ,integral_abs,eps3,eta3, left,F1,F2);
middle_integ:=simpson3point(x0+dx3,dx3,middle_integ, integral_abs,eps3,eta3,F2,middle,F4);
right_integ:=simpson3point(x0+2.0*dx3,dx3,right_integ,integral_abs,eps3,eta3,F4,F5,right);
simpson3poin:=left_integ+middle_integ+right_integ;
end;
Dec(nest_level);
End; {simpson3point}
Begin {adaptive_simpson}
nest_level:=1;
highest_level:=1;
no_evaluations:=3;
adaptive_simpson:=simpson3point(x0,x1-x0,0.0,0.0,eps,eta,F(x0),F(x0+0.5*(x1-x0)),F(x1));
end;{adaptive_simpson}


1.3. Метод Ромберга и его реализация на языке Pascal.

Интегрирование следующим методом – методом Ромберга – основано на правиле трапеций, использующем кусочно-линейное приближение для интегрируемой функции. Основной факт относительно погрешности в методе трапеций следующий.

Теорема. Пусть F(x) – гладкая функция на интервале [a,b], и этот интервал делится на т равных частей, каждая длиной h = . Пусть I(h) обозначает соответствующее приближение метода трапеций:
I(h) = ,
где f_i=F(a+jh) – значение интегрируемой функции в точке a+jh.
Тогда
,
Где a_k – некоторая константа.
Основное здесь то, что погрешность в методе трапеций может быть выражена рядом по четным степеням шага интегрирования h. Построим таблицу значений T_ik:


В нулевой строке T_0k = I((b – a)/2k), так что T₀₀,T₀₁,… являются последовательными приближениями метода трапеций для интеграла, каждое с удвоенным по сравнению с предыдущим числом интервалов. Согласно приведенной выше теореме,
,
где h = ((b – a)/2k.
Отсюда следует, что
,
Поэтому положим
.
В общем случае строим j-ю строку таблицы Ромберга по формуле
,
а оценка погрешности имеет вид
,
где h = (b – a)/2k.
Для работы понадобится не целая таблица, а только последняя вычисленная строка. Число точек выборки на каждом шаге удваивается. Обратите внимание на то, что функцию следует вычислять только в новых точках, которые являются средними точками предыдущих подынтервалов:
F₀ + 2F₁ + 2F₂ + …+ 2F_2n-1 + F_2n =
= (F₀ + 2F₂ + 2F₄ + …+ 2F_2n-2 + F_2n) + 2(F₁ + F₃ +…+F_2n-1).

Таким образом, чтобы модифицировать предыдущее приближение, необходимо вычислить сумму значений функции в новых средних точках. Это делается в цикле со счетчиком k. Метод Ромберга реализован в функции romberg (листинг1.5).