5. Вычисление интегралов с заданной точностью
Уменьшая шаг, мы добиваемся большей точности результата. Это подсказывает алгоритм двойного пересчета, когда вычисления проводятся с шагом h, а затем с шагом h/2. Результаты сравниваются и если разность по модулю полученных квадратур удовлетворяет условию
,
где наперед заданная точность, то расчет заканчивается.
Детальнее проанализировав этот алгоритм, нетрудно заметить, что выполняются лишние действия. Действительно, когда проводится повторный расчет с уменьшенным вдвое шагом, происходит повторное вычисление значений функции в некоторых из узлов. Проиллюстрируем сказанное на методе трапеций.
1. Пусть n=1, т.е. на интервале [a,b] всего одна трапеция. Обозначим квадратурную формулу как
,
где
Если n=2, то
где
Далее при n=4:
,
где
.
Эти выкладки позволяют составить следующую рекуррентную формулу
,
(7.14)
где
;
.
Она позволяет сделать алгоритм двойного
пересчета более экономичным.
Вычисляя квадратурные формулы с шагом h и h/2, можно существенно улучшить результат, составляя комбинацию этих результатов, как это делалось в методах численного дифференцирования. Запишем
;
,
откуда нетрудно
получить постоянную
,
подставляя которую во второе из этих
равенств, получим соотношение
(7.15)
Запишем его
для метода трапеций в узлах
:
подставляя в (7.15 ), получим
формулу Симпсона для частичного интервала. Следовательно, комбинируя в методе двойного пересчета формулы трапеций, получаем метод парабол, имеющий точность на два порядка выше.
Обобщение этих выводов связано с формулой Ромберга:
,
(7.16)
где
,
а результаты удобно представить в виде
треугольной таблицы
где индекс i описывает изменение по строкам, а j - по столбцам.
Процедура заканчивается, когда диагональные элементы удовлетворяют условию:
(7.17)
Составим подпрограмму-функцию, реализующую метод Ромберга
function TrRomb( a,b,eps: real; f1:fun): real;
const m=10;
var q:array{1..m,1..m] of real;
h,s:real; i,j,k,m1,m2: integer;
begin
h:=b-a; i:=1; m1:=1;
q[1,1]:=0.5*h*(f1(a)+f1(b));
repeat
s:=0; Inc(i);
if (i>=m) then
begin writeln(‘Error!’); Exit end;
for k:=1 to m1 do
s:=s+ f1(a+h*(k-0.5));
m1:=m1*2; m2:=1;
q[i,1] := 0.5*(q[i-1,1]+ h*s); h:=0.5*h;
for j:=2 to i do
begin m2:=m2*4;
q[i,j]:=(m2*q[i.j-1] - q[i-1.j-1])/(m2-1)
end;
until (abs(q[i.j] - q[i-1.j-1]) < eps);
TrRomb:= q[i.i]
end;
Приведем
результаты расчетов в виде треугольной
таблицы для выбранного нами теста
0.0000000000
1.5707963268 2.0943951024
1.8961188979 2.0045597550 1.9985707318
1.9742316019 2.0002691699 1.9999831309 2.0000055500
1.9935703438 2.0000165910 1.9999997525 2.0000000163 1.9999999962
1.9983933610 2.0000010334 1.9999999962 2.0000000001 2.0000000000 2.0000000000
Первый столбец – это уже известные нам результаты, когда в методе трапеций мы просто уменьшали шаг вдвое. Видно, что скорость сходимости в таком подходе невелика. Однако смещение по таблице вправо, т.е. комбинация с соответствующим весом предыдущих результатов говорит сама за себя.
Дальнейшее развитие этого алгоритма, связанное с изменением шага интегрирования в зависимости от функциональных колебаний, именуется адаптивной квадратурой. О нем Вы можете прочитать в специальной литературе.