Главная       Скачать       Коммерческая поддержка       FAQ       Forum       О нас       Английская версия

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

Эта страница содержит описание алгоритмов одномерного интегрирования, входящих в состав пакета 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, в отличие от других библиотек, использует для решения этой задачи обратную коммуникацию. Когда требуется вычислить значение функции (или её производных), состояние алгоритма сохраняется в специальной структуре, после чего управление возвращается в вызвавшую программу, которая осуществляет все вычисления и снова вызывает вычислительную подпрограмму.

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

  1. Подготовка структуры данных AutoGKState при помощи одной из подпрограмм инициализации алгоритма (AutoGKSmooth, AutoGKSmoothW, AutoGKSingular).
  2. Вызов подпрограммы AutoGKIteration.
  3. Если подпрограмма вернула False, работа алгоритма завершена и интеграл найден (сам интеграл может быть получен при помощи подпрограммы AutoGKResults).
  4. Если подпрограмма вернула 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 позволит нам избежать ошибки деления на ноль при вычислении значения функции. Пример работы с этими полями можно увидеть в разделе "описание подпрограмм".

Manual entries

C++ autogk subpackage   
C# autogk subpackage   

This article is intended for personal use only.

Скачать ALGLIB

C#

Исходный код на C#

Downloads page

 

C++

Исходный код на C++

Downloads page

 

C++, арифметика высокой точности

Исходный код на C++, использующий библиотеки MPFR/GMP.

Исходный код GMP доступен на сайте gmplib.org. Исходный код MPFR доступен на сайте www.mpfr.org.

Downloads page

 

FreePascal

Исходный код на Free Pascal.

Downloads page

 

Delphi

Исходный код на Delphi.

Downloads page

 

VB.NET

Исходный код на VB.NET.

Downloads page

 

VBA

Исходный код на VBA.

Downloads page

 

Python

Исходный код на Python (CPython и IronPython).

Downloads page

 

 

ALGLIB® - numerical analysis library, 1999-2017.
ALGLIB is registered trademark of the ALGLIB Project.