Участник:Valkyrie/Адаптивное интегрирование

Адаптивное интегрирование

Автор: В.А. Мастинен

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Адаптивное интегрирование - это метод численного интегрирования, в котором интеграл функции [math]f(x)[/math] аппроксимируется с использованием статических квадратурных формул на адаптивно уточненных подынтервалах области интегрирования. Как правило, адаптивные алгоритмы столь же эффективны, как и традиционные алгоритмы для «достаточно гладких» подынтегральных функций, но также неэффективны для «плохих» подынтегральных функций, на которых традиционные алгоритмы терпят неудачу.

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

На практике нередко встречаются случаи, когда подынтегральная функция меняется по-разному на отдельных участках отрезка интегрирования. Это обстоятельство требует такой организации численного алгоритма, при котором он автоматически приспосабливался к характеру изменения функции -вводил разные значения шага интегрирования на отдельных участках отрезка интегрирования. Это позволяет уменьшить машинное время без потери точности результатов расчета. Такой подход используется обычно при задании подынтегральной функции [math]f(x)[/math] в виде формулы, а не в табличном виде.

1.2 Математическое описание алгоритма

Исходные данные : интегрируемая функция [math]f(x)[/math], определенная на отрезке [math][a,b][/math], допустимая погрешность интегрирования [math]\varepsilon[/math].

Вычисляемые данные : приближенное значение [math]Q[/math] интеграла [math]\int\limits_a^b f(x)\,dx[/math].

Пусть [math]Q(a,b)[/math] - квадратурная формула, аппроксимирующая значение интеграла на отрезке [math][a,b][/math].

Оценим [math]\tau \approx |Q(a,b) - \int\limits_a^b f(x)\,dx|[/math].

Если [math]\tau \gt \varepsilon[/math], поделим отрезок [math][a,b][/math] пополам и применим этот же алгоритм к каждому из отрезков [math][a,c][/math] и [math][c,b][/math], где [math] c = \frac{a+b}2 [/math], задав погрешность интегрирования на каждом из отрезков [math]\frac{\varepsilon}2[/math].

Если в качестве квадратурной формулы выбрать формулу Симпсона, то [math]\tau[/math] можно положить равным [math]\frac{|Q(a,c) + Q(c,b) - Q(a,b)|}{15}[/math] ().

Далее в качестве квадратурной формулы будет использоваться формула Симпсона.

1.3 Вычислительное ядро алгоритма

Основную часть алгоритма составляют вычисления квадратурной формулы, следовательно, и значений интегрируемой функций в трех точках, на каждом отрезке интегрирования. Так же на каждом отрезке подсчитывается ошибка интегрирования, которая представляет собой разность обычной квадратурной формулы [math]Q(a,b)[/math] и уточненной [math]\tilde{Q}(a,b) = Q(a,c) + Q(c,b)[/math].

1.4 Макроструктура алгоритма

Как уже записано в описании ядра алгоритма, основную часть алгоритма составляют вычисление квадратурной формулы и подсчет ошибки интегрирования на текущем отрезке интегрирования. После чего решается уточнять ли значение квадратурной формулы на данном отрезке или нет.

1.5 Описание схемы реализации последовательного алгоритма

Последовательный алгоритм представляет собой рекурсивный подсчет приближенного значения интеграла на отрезке с заданной допустимой погрешностью. Начальным отрезком является весь отрезок интегрирования. Далее, на каждом шаге рекурсии происходит подсчет трех квадратурных формул - [math]Q(a,b)[/math], [math]Q(a,c)[/math] и [math]Q(c,b)[/math], после чего оценивается погрешность интегрирования [math]\tau = \frac{|Q(a,c) + Q(c,b) - Q(a,b)|}{15}[/math]. Если [math]\tau \lt \varepsilon[/math], тогда приближенным значением интеграла на данном отрезке считаем величину [math]Q(a,c) + Q(c,b)[/math]. Иначе рекурсивно считаем таким же образом приближенные значения интеграла на отрезках [math][a,c][/math] и [math][c,b][/math], задав допустимой погрешностью на каждом отрезке величину [math]\frac{\varepsilon}{2}[/math], а их сумму примем за результат.

1.6 Последовательная сложность алгоритма

Число операций, необходимых для вычисления приближенного значения интеграла, вообще говоря, зависит от допустимой погрешности интегрирования [math]\varepsilon[/math] и от степени гладкости интегрируемой функции [math]f(x)[/math]. Пусть алгоритм разбил [math]k[/math] отрезков пополам, тогда в окончательном разбиении отрезка интегрирования получилось [math]k+1[/math] отрезков. На каждом таком отрезке и на этапе разбиения приходится считать [math]Q(a,b)[/math], [math]Q(a,c)[/math] и [math]Q(c,b)[/math], то есть необходимо произвести вычисление значений [math]f(x)[/math] в 5 точках, для подсчета каждой квадратуры необходимо выполнить 3 сложения, 1 вычитание, 1 умножение и 2 деления, а также 1 сложение, 1 вычитание и 1 деление для проверки критерия разбиения отрезка пополам и 1 деление при каждом разбиении отрезка. Таким образом, получается:

• [math]5(2k+1)[/math] вычислений значения [math]f(x)[/math] в некоторых точках
• [math]6(2k+1) [/math] сложений и вычитаний
• [math]9k+4[/math] делений и умножений

1.7 Информационный граф

На рисунке 1 приведена работа алгоритма при [math]\varepsilon = 4*10^{-6}[/math], [math]a = 0[/math], [math]b = 1[/math] и [math]f(x) = \left\{ \begin{matrix} sin(\pi x) \text{ при } x \in [0, \frac{1}{2}] \text{,} \\ 1 - 4 (x - \frac{1}{2})^2 \text{ при } x \in [\frac{1}{2}, 1] \text{.}\end{matrix} \right. [/math]

Рисунок 1. Квадратики - задачи интегрирования [math]f(x)[/math] на отрезке [math][a,b][/math] при допустимой погрешности [math]\varepsilon[/math], разбиение задачи на 2 происходит при условии [math]\tau \gt \varepsilon[/math]. Оранжевым выделены операции вычисления квадратурной формулы на данном отрезке интегрирования, желтым - уточненной квадратурной формулы. Зеленым отмечены операции оценивания ошибки интегрирования. Фиолетовым выделены операции сложения двух чисел. Красным отмечен результат работы алгоритма.

Напомним, что уточненная квадратурная формула вычисляется следующим образом - [math]\tilde{Q}(a,b) = Q(a,\frac{a+b}{2}) + Q(\frac{a+b}{2},b)[/math]