Квадратуры Гаусса-Кронрода: адаптивное интегрирование с контролем ошибки
Идеальная формула и проблема оценки погрешности
В области численного интегрирования (квадратур) методы Гаусса являются абсолютным математическим совершенством. В то время как метод Симпсона по 3 равномерным узлам точно интегрирует полиномы только 3-й степени, квадратура Гаусса-Лежандра по тем же 3 узлам способна абсолютно точно проинтегрировать полином 5-й степени! Сдвигая узлы с равномерной сетки в корни ортогональных многочленов Лежандра, метод Гаусса выжимает теоретический максимум алгебраической точности (степень 2n-1 для n узлов).
Однако на практике вычислительной математики у классических формул Гаусса выявился фатальный инженерный недостаток: отсутствие дешевого способа контроля погрешности. Программа-решатель не знает, аналитическую функцию какой сложности ей «скормили». Чтобы гарантировать заданную пользователем точность, алгоритму нужно вычислить интеграл по n узлам, затем удвоить количество узлов (взять 2n) и сравнить результаты. Но узлы Гаусса для разных n не пересекаются (за исключением нулевой точки для нечетных n). Это значит, что вычисленные ранее значения функции придется выбросить в мусорную корзину и вычислять тяжелую подынтегральную функцию во всех новых 2n точках заново. При многократном измельчении это приводит к колоссальным потерям машинного времени.
Гениальное расширение Александра Кронрода
Решение этой сложнейшей проблемы было найдено в 1964 году советским математиком Александром Семеновичем Кронродом. Его идея была настолько изящна, что навсегда вошла в золотой фонд мировой вычислительной математики. Кронрод поставил задачу: как добавить к уже существующим n узлам Гаусса оптимальное количество новых узлов, чтобы старые вычисления (узлы) были полностью переиспользованы, а точность новой, составной формулы была бы максимально возможной?
Кронрод математически доказал, что оптимальным является добавление строго n+1 новых узлов (они лежат между узлами Гаусса). Полученная составная квадратура по 2n+1 узлам обладает алгебраической степенью точности 3n+1. Теперь алгоритм работает так: сначала вычисляется грубый интеграл по n узлам Гаусса. Затем вычисляются значения функции только в n+1 новых точках Кронрода. Старые и новые значения складываются с новыми весами, давая высокоточный интеграл Кронрода. Разница между интегралом Гаусса и интегралом Кронрода дает великолепную, практически бесплатную оценку текущей ошибки интегрирования!
Библиотека QUADPACK и адаптивное разбиение
Самым популярным стандартом в индустрии стала пара узлов: 7-точечная формула Гаусса и 15-точечная формула Кронрода (G7/K15). На базе этой элегантной пары бельгийскими математиками (Писсенс и де Донкер) был написан пакет программ QUADPACK.
Этот пакет реализует глобальную адаптивную стратегию. Алгоритм вычисляет интеграл на всем заданном отрезке методом Гаусса-Кронрода и оценивает ошибку. Если ошибка превышает заданный допуск, алгоритм не увеличивает порядок полиномов, а разрезает отрезок ровно пополам (строит бинарное дерево). В каждой половине процедура повторяется. Алгоритм поддерживает список (очередь) всех подотрезков, сортируя их по величине вносимой ошибки, и агрессивно дробит только те участки оси, где функция ведет себя нерегулярно (резкие скачки, осцилляции). Этот мощный, надежный и экономный алгоритм (qag в QUADPACK) сегодня является скрытым мотором интегрирования в функциях `integral` системы MATLAB и `scipy.integrate.quad` в экосистеме Python.