Квадратуры Гаусса-Лежандра очень точны, но НАСКОЛЬКО они точны? Разумеется, есть формулы, по котором можно получить верхнюю границу для погрешности, но на практике их не используют. Скажем, для формулы с двадцатью шестью узлами надо оценить сверху пятьдесят вторую (!) производную от интегрируемой функции, чтобы получить оценку погрешности.

Можно попробовать следующий способ, дающий неплохие результаты: вычислить интеграл по формуле с M узлами, а потом вычислить второй, по формуле с 2M узлами, и сравнить их. Если функция достаточно гладкая, то второй интеграл будет намного точнее первого, и их разность будет близка к погрешности первого интеграла. Ну а погрешность второго интеграла будет ещё меньше.

Есть только один недостаток - слишком уж много раз требуется вычислять интегрируемую функцию. Поэтому гораздо эффективнее использовать квадратурную формулу Гаусса-Кронрода.

Квадратурная формула Гаусса-Кронрода

Идея этой формулы проста- пусть мы вычислили первый интеграл, с M узлами, по квадратурной формуле Гаусса-Лежандра. Если мы станем вычислять интеграл с 2M узлами по этой же формуле, то новые узлы интегрирования не будут совпадать со старыми, и ранее вычисленные значения функции пропадут, не смогут быть использованы.

Вместе с тем, есть способ построения формулы с 2M+1 узлами, которая получается путем добавления к старой формуле с M узлами, (M+1)-ой новой точки. Т.е. при интегрировании по более точной формуле мы можем использовать значения, полученные при интегрировании по менее точной формуле. К сожалению, формула Гаусса-Кронрода с 2M+1 узлами менее точна, чем формула Гаусса-Лежандра с тем же количеством узлов: порядок первой - 3M+1, порядок второй - 4M+1. И всё же дело того стоит.

Алгоритм вычисления узлов и весов

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

Погрешность интегрирования

К сожалению, у меня отсутствует информация о погрешности данной формулы. Если она есть у вас, буду рад, если вы мне её пришлете.

Описание модуля

В модуле содержатся две подпрограммы: BuildGaussKronrodQuadrature и IntegralGaussKronrod. Первая по заданному N выводит из таблицы набор узлов и весов формулы Гаусса-Кронрода, вторая позволяет вычислить интеграл от функции, используя результаты работы первой.

Подпрограмма BuildGaussKronrodQuadrature

На входе алгоритм получает:

• N - требуемое число точек квадратуры Кронрода (т.е. второй квадратуры, с N = 2M+1 узлами). Допустимые значения N: 15, 21, 31, 41, 51, 61.

На выходе алгоритм возвращает:

• X - массив узлов формул Гаусса (M узлов) и Кронрода (2M+1 узел). Нумерация элементов от 0 до N-1. Элементы с нечетными номерами принадлежат и формуле Гаусса, и формуле Кронрода. Элементы массива с четными номерами принадлежат только формуле Кронрода.
• WKronrod - массив весовых коэффициентов формулы Кронрода. Нумерация элементов от 0 до N-1.
• WGauss - массив весовых коэффициентов формулы Гаусса. Нумерация элементов от 0 до N-1. Элементы массива WGauss с четными номерами равны 0 и в формулу Гаусса не входят.

Подпрограмма IntegralGaussKronrod

На входе алгоритм получает:

• A, B - отрезок интегрирования
• N - число точек в квадратуре Гаусса-Кронрода. Допустимые значения N: см. описание BuildGaussKronrodQuadrature.

На выходе алгоритм возвращает:

• Err - оценка сверху абсолютной погрешности интегрирования. Нет никаких гарантий, что погрешность окажется меньше оценки, но для гладких функций погрешность обычно на несколько порядков меньше своей оценки.

Результат: интеграл F на отрезке [A, B]

Report a bug