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

Эта страница содержит описание алгоритмов одномерного интегрирования, входящих в состав пакета ALGLIB.

Адаптивная квадратура: интегрирование гладкой функции

В подпрограмме AutoGKSmooth реализован простейший алгоритм адаптивного интегрирования, позволяющий находить интегралы гладких функций. В основе алгоритма лежит рекурсивное разбиение интервала интегрирования и применение на получившихся субинтервалах квадратурной формулы Гаусса-Кронрода (в текущей реализации используется формула с 15 узлами). Эта квадратурная формула позволяет вычислить одновременно и сам интеграл, и оценку его погрешности. Субинтервалы объединяются в кучу (heap), после чего из кучи извлекается интервал с наибольшей абсолютной ошибкой и снова подвергается разбиению. Процесс повторяется, пока суммарная относительная погрешность не приблизится к машинной точности (текущая реализация использует следующий критерий - |Err| ≤ C·ε·|Integral|).

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

Адаптивная квадратура: функции с интегрируемыми сингулярностями на концах интервала

Рассмотренные выше подпрограммы хорошо справляются с гладкими функциями. Их эффективность заметно снижается, если функция имеет интегрируемые сингулярности на одном из концов отрезка, т.е. ведет себя, как (x-a)^α или (b-x)^β (сингулярность интегрируема, если α > -1 и β > -1). Иногда требуются тысячи или десятки тысяч разбиений, чтобы погрешность интегрирования стала достаточно мала.

Подпрограмма AutoGKSingular решает эту проблему. Используя замену переменных, она может сделать легко интегрируемой любую сингулярность вида (x-a)^α или (b-x)^β - при условии, что известны коэффициенты α и β. Впрочем, если их точная величины неизвестны, можно воспользоваться оценкой снизу - с условием, что оценки коэффициентов тоже строго больше, чем -1.

Обратная коммуникация и вычисление значений функции

Алгоритм интегрирования в ходе своей работы должен получать значения функции в выбранных им точках. В большинстве программных пакетов эта проблема решается путем передачи указателя на функцию (C++, Delphi) или делегата (C#), который осуществляет эту операцию.

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

Таким образом, работа с алгоритмом оптимизации осуществляется в следующей последовательности:

Подготовка структуры данных AutoGKState при помощи одной из подпрограмм инициализации алгоритма (AutoGKSmooth, AutoGKSmoothW, AutoGKSingular).
Вызов подпрограммы AutoGKIteration.
Если подпрограмма вернула False, работа алгоритма завершена и интеграл найден (сам интеграл может быть получен при помощи подпрограммы AutoGKResults).
Если подпрограмма вернула True, подпрограмма требует информацию о функции. Вычислите значение функции в точке State.X и сохраните его в поле State.F, после чего повторно вызовите подпрограмму AutoGKIteration.

Следует отметить, что при использовании подпрограммы AutoGKSingular (т.е. при интегрировании функций с сингулярностями на концах интервала) может потребоваться значение функции в точке, расположенной слишком близко к одному из концов отрезка. Например, интегрируя функцию f(x) = (1-x)^-0.95 на отрезке [0, 1], нам может потребоваться значение f(1-10^-25). Очевидно, что при использовании арифметики обычной точности произойдет потеря точности, 1-10^-25 будет округлено до 1, что приведет к ошибке деления на ноль.

Для того, чтобы избежать такой ситуации, предприняты следующие меры. Когда требуется вычислить значение функции, алгоритм заполняет не только поле State.X, но также и поля State.XMinusA и State.BMinusX. Если рассмотреть пример, приведенный выше, то поля State.X и State.XMinusA будут равны 1 (при использовании арифметики двойной точности), поле State.BMinusX будет равно 10^-25. Использование поля State.BMinusX позволит нам избежать ошибки деления на ноль при вычислении значения функции. Пример работы с этими полями можно увидеть в разделе "описание подпрограмм".