Прецизионный расчет одномерных квадратур

Аннотация

Вычисление квадратур возникает во многих физических и технических приложениях. Отобраны лучшие квадратурные формулы и проведено их количественное сравнение на ряде представительных примеров.

Полный текст

ВВЕДЕНИЕ

Во многих прикладных задачах требуется экономичный расчет квадратур

$I=\underset{0}{\overset{1}{\int }}f\left(x\right)dx$ (1)

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

В литературе описано большое количество методов вычисления таких квадратур. Наиболее употребительны формулы трапеций, средних и Симпсона на равномерных или неравномерных сетках, квадратуры Гаусса–Кристоффеля и гауссово-сеточный метод [1]. Для улучшения сходимости применяют различные замены переменных интегрирования [2], а также используют поправки Эйлера–Маклорена [1, 3].

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

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

ВЫЧИСЛЕНИЕ КВАДРАТУР

Квадратуры трапеций, средних и Симпсона имеют степенную скорость сходимости. Погрешность зависит от шага h сетки как O(h^m). Экстраполяционное уточнение по Ричардсону и уточнение по формулам Эйлера–Маклорена повышают порядок точности m, однако закон сходимости остается степенным. Эти методы общеизвестны, поэтому в данной работе они не рассматриваются.

Кардинально точнее методы, имеющие сверхстепенную либо экспоненциальную сходимость. Их погрешность зависит от шага сетки по закону, близкому к ∼ exp(–1 / h). Это означает, что при уменьшении шага вдвое число верных знаков в ответе примерно удваивается. Такую скорость сходимости реализуют следующие методы.

1. Квадратура Гаусса–Кристоффеля.
2. Сеточные формулы средних и трапеций для периодической подынтегральной функции. повысить, используя процедуру экстраполяции точности. Соответствующие квадратуры называются экстраполяционными.
3. Сеточные формулы средних и трапеций для непериодической подынтегральной функции после специальной замены переменных.

КВАДРАТУРЫ ГАУССА–КРИСТОФФЕЛЯ

Эти квадратуры хорошо известны. Напомним, что они строятся для интегралов вида

$I=\underset{0}{\overset{1}{\int }}u\left(x\right)\rho \left(x\right)dx$ (2)

с весовой функцией ρ(x). Квадратура имеет вид

$I_N=\sum _{j=1}^N{c}_{j}u\left(x_j\right)$. (3)

Здесь x_j и c_j – узлы и веса, соответственно. Тем самым, квадратура (3) содержит 2N параметров. Они выбираются так, чтобы квадратура была точна для полинома степени 2N – 1.

Квадратуры Гаусса–Кристоффеля сталкиваются с двумя трудностями. Во-первых, нахождение узлов и весов этих квадратур представляет отдельную проблему. Эти узлы и веса зависят от вида ρ(x) и найдены только для отдельных частных случаев.

Во-вторых, в оценку погрешности квадратуры, содержащей N узлов, входит производная u^(2N)(x). Поэтому для реализации экспоненциальной сходимости необходимо, чтобы все производные были ограничены на отрезке интегрирования. Если u(x) имеет особенность (например, неограниченную производную), то эту особенность выделяют в весовую функцию ρ(x). Это вынуждает строить узлы и веса для каждой задачи в отдельности (см., например, [4]).

Мы реализовали алгоритм [5] вычисления узлов и весов квадратуры Гаусса–Кристоффеля для единичной весовой функции ρ(x) = 1. Это пакет программ на языке среды Matlab/Octave, выложенный в открытый доступ [6]. Его отличительной особенностью является то, что узлы и веса вычисляются с произвольной наперед заданной разрядностью чисел. Это является преимуществом данного пакета по сравнению с аналогичными (см., например, [5]).

На практике также употребителен гауссово-сеточный метод. Он состоит в том, что отрезок интегрирования разбивается на шаги длины h, для каждого из которых записывается квадратура Гаусса с небольшим числом узлов N = 2 – 5. Этот метод предъявляет такие же требования к гладкости функции, как классическая квадратура Гаусса.

ЭКСТРАПОЛЯЦИОННАЯ КВАДРАТУРА

Пусть функция f(x) периодична с периодом 1. Заменой переменных z = exp(2πix) перейдем к интегралу по единичной окружности

$I=\frac{1}{2\pi i}\underset{\left|z\right|=1}{\int }g\left(z\right)dz, g\left(z\right)=\frac{1}{z}f\left(\ln z\right)$, (4)

и запишем для него формулу трапеций

$I_N=\frac{1}{N}\sum _{n=1}^{N}g\left(z_n\right)z_n$. (5)

Для этой квадратуры в [7, 8] построены асимптотически точные оценки погрешности. Скорость сходимости определяется положением особых точек функции g(z) на комплексной плоскости.

В качестве примера рассмотрим простейший случай, когда g(z) имеет один полюс кратности q в точке z = a внутри единичной окружности. Тогда g(z) представима в виде

$g\left(z\right)=\frac{h\left(z\right)}{{(z-a)}^q}$, (6)

где числитель h(z) не имеет особых точек, кроме, может быть, бесконечно удаленной. В этом случае оценка погрешности квадратуры (5) имеет вид

${\delta }_N=\frac{1}{(q-1)!}\frac{{\partial }^{q-1}}{\partial a^{q-1}}\left[h\left(a\right)\left(1-\frac{1}{1-a^N}\right)\right]$. (7)

В [8] предложен ряд обобщений оценки (7). Функция g(z) может иметь произвольное конечное число полюсов произвольного порядка. Положения и кратности полюсов могут быть неизвестны заранее. В этом случае они вычисляются в ходе сгущения сеток. Оценки, аналогичные (7), построены для квадратуры средних.

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

${\tilde{I}}_N=I_N+{\delta }_N$. (8)

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

ЗАМЕНА ПЕРЕМЕННЫХ

Пусть теперь f(x) не является периодической. Запишем для интеграла (1) формулу средних

$I_N=\frac{1}{N}\sum _{n=1}^{N}f\left(x_{n-1/2}\right)$. (9)

Выполним замену переменной [9]

$t\left(\xi \right)=\frac{(\xi -0.5)}{\xi (1-\xi )}, x\left(t\right)=\frac{1}{2}+\frac{1}{2}th t$. (10)

Тогда интеграл примет вид

$I=\underset{0}{\overset{1}{\int }}f\left(x\right(\xi \left)\right)x_{\xi }\left(\xi \right)d\xi \equiv \underset{0}{\overset{1}{\int }}\tilde{f}\left(\xi \right)d\xi$. (11)

Легко видеть, что множитель x_ξ быстро стремится к нулю при ξ → 0 + 0 и ξ → 1 – 0. Также нетрудно убедиться, что $\tilde{f}\left(\xi \right)$ стремится к нулю вместе со всеми производными при стремлении к точкам x = 0,1 изнутри отрезка. Поэтому $\tilde{f}\left(\xi \right)$ допускает бесконечно гладкое периодическое продолжение за пределы отрезка интегрирования.

Согласно формуле Эйлера–Маклорена, степенная часть погрешности квадратуры средних и трапеций содержит слагаемые вида $N^{-2m}\left({\tilde{f}}^{(2m-1)}\right(1)-{\tilde{f}}^{(2m-1)}(0\left)\right)$. Эти слагаемые обращаются в нуль, и сходимость становится сверхстепенной.

Отметим, что такой закон сходимости реализуется независимо от поведения f(x) вблизи концов отрезка (0,1): она может иметь неограниченную производную и даже обращаться в бесконечность (несобственный интеграл).

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

СРАВНЕНИЕ КВАДРАТУР

Проведем расчеты тестовых примеров и сравним описанные выше квадратурные формулы. Рассмотрим сначала периодическую подынтегральную функции. Пусть

$I=\frac{1}{2\pi i}\underset{\left|z\right|=1}{\int }\frac{1}{(z-a_1)(z-a_2)}th\frac{1}{z-b}dz$. (12)

Положим a₁ = 2 / π, a₂ = π / 2, b = 1000 / π. Точное значение этого интеграла нетрудно найти с помощью вычетов.

На рис. 1 показаны графики погрешности для разных квадратурных формул. Масштаб полулогарифмический, так что экспоненциальной сходимости соответствует прямая. Кружками обозначена погрешность формулы трапеций. Треугольники соответствуют погрешности формулы Гаусса–Кристоффеля. Темные квадраты – погрешность экстраполяционной квадратуры, в которой уточнение проводится при известных положениях полюсов. Светлые квадраты – погрешность экстраполяционной квадратуры, в которой уточнение проводилось при заранее неизвестных положениях полюсов.

Рис. 1. Погрешности расчета тестового интеграла (12), обозначения см. текст

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

Рассмотрим второй пример, в котором подынтегральная функция является непериодической. Пусть

$I=\underset{0}{\overset{1}{\int }}{x}^{\beta -1}{e}^{-x}dx, \beta =1.5$. (13)

Точный ответ есть нижняя неполная гамма-функция γ(1, β). Подынтегральная функция имеет ограниченную гладкость: производная f' ∼ x^–0.5 неограниченно возрастает вблизи x = 0.

На рис. 2 показаны графики сходимости разных квадратур в полулогарифмическом масштабе. Здесь светлыми кружками обозначена погрешность обыкновенной формулы средних. Темные кружки соответствуют погрешности формулы средних с заменой переменных (10). Треугольники – погрешность квадратуры Гаусса–Кристоффеля.

Рис. 2. Погрешности расчета тестового интеграла (13), обозначения см. текст

Видно, что замена переменных (10) кардинально ускоряет сходимость. Формула средних с заменой переменных значительно превосходит по точности формулу Гаусса–Кристоффеля. Выигрыш достигает 8 порядков.

ЗАКЛЮЧЕНИЕ

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

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

Если подынтегральная функция периодична, но ближайшая особая точка является существенно особой либо точкой ветвления, то наиболее точной оказывается квадратура Гаусса–Кристоффеля с единичным весом. Для ее построения целесообразно использовать пакет программ [6].

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

Во всех остальных случаях целесообразно использовать формулу средних с заменой переменных (10).

Вторая точка зрения заключается в том, что предпочтительнее использовать единообразную стратегию расчета для как можно более широкого круга задач. В этом случае мы рекомендуем формулу средних с заменой переменных (10). Если подынтегральная функция имеет особенности, то исходный отрезок интегрирования необходимо разбить на промежуточные отрезки, на которых подынтегральная функция является бесконечно гладкой. Далее на каждом таком подотрезке нужно ввести замену переменных, аналогичную (10), и записать формулу средних с равномерным шагом.

Работа выполнена при поддержке Российского научного фонда (проект № 22-71-00028, https://rscf.ru/project/22-71-00028).