Известия РАН. Серия физическая, 2022, T. 86, № 11, стр. 1628-1633
Квадратуры со сверхстепенной сходимостью
М. А. Тинтул 1, *, В. С. Хохлачев 1, А. А. Белов 1, 2
1 Федеральное государственное бюджетное образовательное учреждение высшего образования
“Московский государственный университет имени М.В. Ломоносова”, Физический факультет
Москва, Россия
2 Федеральное государственное автономное образовательное учреждение высшего образования
“Российский университет дружбы народов”
Москва, Россия
* E-mail: maksim.tintul@mail.ru
Поступила в редакцию 30.06.2022
После доработки 15.07.2022
Принята к публикации 22.07.2022
- EDN: HDNZDQ
- DOI: 10.31857/S0367676522110308
Аннотация
Вычисление квадратур возникает во многих физических и технических приложениях. Предложена замена переменных интегрирования, кардинально повышающая точность формулы средних. Для бесконечно гладких подынтегральных функций закон сходимости становится сверхстепенным. Он существенно быстрее степенного и близок к экспоненциальному. Для подынтегральных функций с ограниченной гладкостью реализуется степенная сходимость с максимально достижимым порядком точности.
ВВЕДЕНИЕ
Во многих физических задачах требуется приближенно вычислять интегралы, не берущиеся в элементарных функциях. Приведем некоторые примеры:
1) Вычисление специальных функций математической физики: функции Ферми–Дирака, равные моментам фермиевского распределения, гамма-функция, цилиндрические функции и ряд других.
2) Расчет фурье-коэффициентов заданной функции, преобразования Фурье и Лапласа.
3) Численное решение интегральных уравнений, как корректно поставленных, так и некорректных.
4) Решение краевых задач для уравнений в частных производных (включая задачи на собственные значения), записанных в интегральной форме и т.д.
Такие интегралы необходимо вычислять с высокой точностью вплоть до ошибок компьютерного округления.
Как правило, для сеточного вычисления квадратур используют методы трапеций, средних и Симпсона на равномерной сетке. Для этих методов хорошо известна мажорантная оценка погрешности. Для формул трапеций и средних она составляет $O({{h}^{2}}),$ для формулы Симпсона – $O({{h}^{4}}),$ где $h$ – шаг сетки. Существуют способы повышения точности: расчет на наборе сгущающихся сеток и экстраполяционное уточнение по методу Ричардсона, уточнение по формуле Эйлера–Маклорена и др. [1, 2]. Все эти способы дают степенную зависимость погрешности от шага сетки $O({{h}^{m}}).$
Если подынтегральная функция является периодической, и интеграл вычисляется по полному периоду, то зависимость погрешности от шага становится не степенной, а экспоненциальной $\sim {\kern 1pt} \exp ( - {1 \mathord{\left/ {\vphantom {1 h}} \right. \kern-0em} h})$ [3–5]. Это означает, что при уменьшении шага вдвое число верных знаков в ответе примерно удваивается. Такая скорость сходимости намного быстрее степенной. Однако соответствующий класс подынтегральных функций довольно узок. В литературе предпринимались попытки расширить этот класс [6–9], однако они были признаны [7] неудачными.
В данной работе предложен подход, кардинально ускоряющий сходимость формулы средних. Он основан на специальной замене переменных интегрирования. Подынтегральная функция может быть непериодической. Если она бесконечно гладкая, то предлагаемая замена обеспечивает сверхстепенную сходимость квадратуры. Такая скорость сходимости существенно быстрее степенной и близка к экспоненциальной.
Если подынтегральная функция имеет ограниченную гладкость, то предлагаемый метод дает степенную сходимость с максимально достижимым порядком точности.
Предлагаемый подход не требует априорной информации о характере подынтегральной функции и единообразно применим к широкому кругу задач. Класс подынтегральных функций, для которых реализуется сверхстепенная сходимость квадратур, существенно расширяется.
ЗАМЕНА ПЕРЕМЕННЫХ ИНТЕГРИРОВАНИЯ
Рассмотрим интеграл
Выполним замену переменных интегрирования в два этапа. Сначала с помощью дробно-полиномиального преобразования $t(x)$ отобразим отрезок $x \in (0,1)$ на прямую $t \in ( - \infty , + \infty ).$ Затем эту прямую отобразим на отрезок $\xi \in (0,1)$ с помощью преобразования $t(\xi ),$ производные которого стремятся к нулю вблизи $\xi = 0$ и $\xi = 1$ быстрее любой степени ${{\xi }^{m}}.$Такие замены можно составить различными способами. В данной работе рассматривалось преобразование
(2)
$t(\xi ) = \frac{{A(\xi - 0.5)}}{{{{\xi }^{\alpha }}{{{(1 - \xi )}}^{\alpha }}}},\,\,\,\,x(t) = \frac{1}{2} + \frac{1}{2}{\text{th(}}Bt{\text{),}}$
После замены (2) интеграл преобразуется к виду
(3)
$I = \int\limits_0^1 {\tilde {f}(\xi )d\xi } ,\,\,\,\,\tilde {f}(\xi ) = f\{ x[t(\xi )]\} {{x}_{t}}[t(\xi )]{{t}_{\xi }}(\xi ).$Периодическое продолжение
Покажем, что новая подынтегральная функция $\tilde {f}(\xi )$ допускает бесконечно гладкое периодическое продолжение за границы отрезка $\xi \in (0,1).$
Выражение ${{t}_{\xi }}(\xi )\sim {{\xi }^{{ - \alpha - 1}}}{{(1 - \xi )}^{{ - \alpha - 1}}}$ имеет полюсы на концах отрезка $\xi = 0$ и $\xi = 1.$ Однако при $\xi \to 0 + 0$ и $\xi \to 1 - 0$ производная ${{x}_{t}}\sim \exp ( - {{\xi }^{{ - \alpha }}}{{(1 - \xi )}^{{ - \alpha }}})$ стремится к нулю существенно быстрее. В результате ${{x}_{t}}{{t}_{\xi }} \to 0$ при стремлении к точкам $\xi = 0$ и $\xi = 1$ изнутри отрезка. Поэтому $\tilde {f}(\xi )$ обращается в нуль на границах отрезка.
Аналогично можно показать, что все производные этой функции стремятся к нулю при $\xi \to 0 + 0$ и $\xi \to 1 - 0.$ Например, первая производная имеет вид
(4)
${{\tilde {f}}_{\xi }} = {{f}_{x}}x_{t}^{2}t_{\xi }^{2} + f{{x}_{{tt}}}t_{\xi }^{2} + f{{x}_{t}}{{t}_{{\xi \xi }}}.$Все производные ${{d{{t}^{m}}} \mathord{\left/ {\vphantom {{d{{t}^{m}}} {d{{\xi }^{m}}}}} \right. \kern-0em} {d{{\xi }^{m}}}}\sim {{\xi }^{{ - \alpha - m}}}{{(1 - \xi )}^{{ - \alpha - m}}}$ на границах отрезка имеют полюсы, которые умножаются на выражение $\sim {\kern 1pt} \exp ( - {{\xi }^{{ - \alpha }}}{{(1 - \xi )}^{{ - \alpha }}})$ в различных степенях. Поэтому при $\xi \to 0 + 0$ и $\xi \to 1 - 0$ имеем ${{\tilde {f}}_{\xi }} \to 0.$ То же верно и для более высоких производных ${{d{{{\tilde {f}}}_{\xi }}^{m}} \mathord{\left/ {\vphantom {{d{{{\tilde {f}}}_{\xi }}^{m}} {d{{\xi }^{m}}}}} \right. \kern-0em} {d{{\xi }^{m}}}}.$ Таким образом, подынтегральная функция $\tilde {f}$ может быть бесконечно гладко периодически продолжена за границы отрезка $\xi \in (0,1).$
СХОДИМОСТЬ ФОРМУЛЫ СРЕДНИХ
На отрезке $\xi \in (0,1)$ введем равномерную сетку с шагом $h = {1 \mathord{\left/ {\vphantom {1 N}} \right. \kern-0em} N}.$ Полуцелые узлы обозначим через ${{\xi }_{{n + 1{\text{/}}2}}} = (n - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2})h,$ $n = 1,...,N.$ Запишем квадратуру средних
Справедлива
Теорема
А) Если $f(x)$ является бесконечно гладкой на отрезке $x \in (0,1),$ то квадратура (5) имеет сверхстепенную сходимость. Б) Если $f(x)$ имеет $j$ непрерывных производных на $x \in (0,1),$ $(j + 1)$-я производная имеет разрыв в точке $x = a \in (0,1),$ и эта точка является узлом сетки, то квадратура (5) имеет степенную сходимость. Порядок точности равен $j + 2,$ если $j$ четно, и $j + 3,$ если $j$ нечетно. Такой порядок точности является максимальным при данной гладкости подынтегральной функции.
Доказательство
Докажем утверждение А). Степенная часть погрешности квадратуры средних (5) описывается формулой Эйлера–Маклорена [1]. Она содержит разности нечетных производных на концах отрезка интегрирования
(6)
$\delta = \sum\limits_{k = 1}^\infty {{{b}_{k}}{{h}^{{2k}}}\left( {{{{\tilde {f}}}^{{(2k - 1)}}}(1) - {{{\tilde {f}}}^{{(2k - 1)}}}(0)} \right)} ,\,\,\,\,{{b}_{k}} = {\text{const}}{\text{.}}$Как отмечалось выше, после замены переменных (2) производные ${{\tilde {f}}^{{(k)}}}(\xi ) \to 0$ при $\xi \to 0 + 0$ и $\xi \to 1 - 0.$ Все слагаемые в сумме (6) обращаются в нуль. Поэтому в погрешности формулы средних не остается степенных членов, и сходимость оказывается сверхстепенной.
Докажем утверждение Б). При указанных предположениях степенной вклад в погрешность формулы средних имеет вид
(7)
$\begin{gathered} \delta = \sum\limits_{k = 1}^\infty {{{b}_{k}}{{h}^{{2k}}}\left( {{{{\tilde {f}}}^{{(2k - 1)}}}(1) - {{{\tilde {f}}}^{{(2k - 1)}}}(0)} \right)} + \\ + \,\,\sum\limits_{k = 1}^K {{{b}_{k}}{{h}^{{2k}}}\left( {{{{\tilde {f}}}^{{(2k - 1)}}}(a - 0) - {{{\tilde {f}}}^{{(2k - 1)}}}(a + 0)} \right)} . \\ \end{gathered} $Первая сумма в (7) аналогична (6). После замены переменных (2) она обращается в нуль.
Вторая сумма есть погрешность, возникающая из-за особенности в точке $a.$ Если $2k - 1 \leqslant j,$ то в силу непрерывности правые и левые предельные значения производных порядка $2k - 1$ одинаковы ${{\tilde {f}}^{{(2k - 1)}}}(a - 0) = {{\tilde {f}}^{{(2k - 1)}}}(a + 0).$ Каков предел суммирования $K?$ Поскольку ${{\tilde {f}}^{{(j + 1)}}}$ разрывна в точке $a,$ и в (7) входят только нечетные производные, то возможны два случая. Если $j$ нечетное, то $2K - 1 = j + 2.$ Тогда $\delta = O({{h}^{{j + 3}}}).$ Если $j$ четное, то $2K - 1 = j + 1,$ и $\delta = O({{h}^{{j + 2}}}).$ Очевидно, такой порядок точности является максимальным, т.е. он не может быть повышен. Теорема доказана.
Замечание
В литературе описаны [6–9] замены переменных, аналогичные (2). В этих работах использовались формулы трапеций и Симпсона, в которых необходимо вычислять подынтегральную функцию в граничных точках. Однако после замены (2) подынтегральная функция $\tilde {f}(\xi )$ имеет существенно особые точки в границах отрезка $\xi = 0$ и $\xi = 1.$ Поэтому вычисление $\tilde {f}(0)$ и $\tilde {f}(1)$ представляет проблему; в частности, возникает переполнение компьютерных чисел.
Чтобы этого избежать, в [6] было предложено обрезать отрезок интегрирования, т.е. вместо $\xi \in (0,1)$ рассматривать $\xi \in (\varepsilon ,1 - \varepsilon ),$ где $\varepsilon $ – некоторое малое число. Такое обрезание внесло существенную погрешность, и реализовать сверхстепенную сходимость не удалось. Авторы работы [7] провели численные эксперименты и установили, что по количественной точности этот подход уступает формуле Симпсона без замены переменных. Поэтому этот подход был признан неперспективным [7].
Мы используем формулу средних, в которой не требуется вычислять $\tilde {f}(0)$ и $\tilde {f}(1).$ Поэтому описанная трудность не возникает, и реализуется сверхстепенная сходимость.
АПРОБАЦИЯ МЕТОДА
Бесконечно гладкая подынтегральная функция
В качестве примера рассмотрим тестовый интеграл с известным точным значением
(8)
$I = \int\limits_0^1 {{{{{e}^{x}}} \mathord{\left/ {\vphantom {{{{e}^{x}}} {(e - 1)}}} \right. \kern-0em} {(e - 1)}}dx} = 1.$Подынтегральная функция является бесконечно гладкой.
Расчет проводился на наборе сеток с различными $N = 2,4,8,...$ На каждой сетке вычислялась квадратура средних и ее погрешность $\Delta = \left| {I - {{I}_{N}}} \right|,$ равная разности численного и точного интегралов. На рис. 2 приведен график погрешности $\Delta $ в зависимости от числа шагов сетки $N.$ Масштаб графика полулогарифмический. В таком масштабе экспоненциальной сходимости соответствует прямая линия, а степенной – логарифмическая кривая.
Темные кружки соответствуют расчету с заменой переменных (2), светлые кружки – расчету без нее. Видно, что предлагаемая замена переменных кардинально повышает точность: уже при $N\sim 100$ погрешность составляет $\Delta \sim {{10}^{{ - 14}}},$ что сопоставимо с ошибками округления. Выигрыш по точности по сравнению с расчетом без замены переменных достигает 10 порядков. Скорость сходимости несколько уступает экспоненциальной, но кардинально превосходит степенную.
Из-за наличия существенно особых точек $\tilde {f}$ на границах отрезка зависимость погрешности от числа шагов является немонотонной и знакопеременной [10, 11]. На данном графике это видно по немонотонному поведению кривой. Локальные минимумы соответствуют смене знака погрешности.
Таким образом, предложенная замена кардинально повышает точность формулы средних. Мы рекомендуем ее к широкому применению в практических вычислениях.
Подынтегральная функция с ограниченной гладкостью
Нередко в приложениях нужно вычислять интегралы от кусочно-заданных сплайн-аппроксимаций и интерполянт. Они имеют ограниченную гладкость. Так, простейшая линейная интерполяция непрерывна, но имеет разрывы первой производной. Кубический сплайн непрерывен вместе со второй производной, а третья производная испытывает разрыв.
В качестве примера рассмотрим интеграл от функции
(9)
$f(x) = \left\{ \begin{gathered} 1,\,\,\,\,x < 0.5, \hfill \\ 1 + {{(2x - 1)}^{m}},\,\,\,\,x \geqslant 0.5, \hfill \\ \end{gathered} \right.$Таблица 1.
Тест (9)
| $m$ | $I$ | $q$ |
|---|---|---|
| $1$ | $1 + 2{{e}^{{0.5}}} - e$ | $2$ |
| $2$ | $1 - 8{{e}^{{0.5}}} + 5e$ | $4$ |
| $3$ | $1 + 48{{e}^{{0.5}}} - 29e$ | $4$ |
| $4$ | $1 - 384{{e}^{{0.5}}} + 233e$ | $6$ |
| $5$ | $1 + 3840{{e}^{{0.5}}} - 2329e$ | $6$ |
Расчет проводился на нескольких сгущающихся сетках. Они выбирались так, чтобы особенность $x = 0.5$ являлась узлом. Например, для этого достаточно брать только четные $N.$ Полученные погрешности в зависимости от числа шагов проведены на рис. 3. Масштаб графика двойной логарифмический. Поэтому степенной сходимости соответствует прямая линия, наклон которой равен порядку точности. Цифры около линий – значения $m.$
Видно, что на достаточно подробных сетках кривые для каждого $m$ стремятся к прямым линиям, т.е. реализуется степенная сходимость. Соответствующие порядки точности $q$ приведены в табл. 1 . Они полностью согласуются с теоремой 1.
На грубых сетках поведение кривых является нерегулярным. Погрешность зависит от $N$ немонотонно, меняет знак и убывает существенно быстрее степенного закона. По-видимому, закон сходимости на грубых сетках соответствует сверхстепенному (аналогично рис. 2). Это наблюдение требует теоретического осмысления.
Для сравнения на рис. 3 приведена погрешность расчета по формуле средних без замены переменных. Для всех рассмотренных $m$ погрешности были примерно одинаковы, поэтому мы показали их одной линией. Она обозначена звездочкой (*). Эта линия соответствует степенной сходимости со вторым порядком точности. Видно, что при $m \geqslant 2$ (т.е. при наличии хотя бы одной непрерывной производной) предложенная замена повышает порядок точности и резко уменьшает количественную погрешность. На рис. 3 выигрыш по точности составил от 3 до 8 порядков.
Мы проводили аналогичный расчет, используя сетки, в которых особенность $x = 0.5$ не попадала в узел. Для этого достаточно брать только нечетные $N.$ Этот случай не подпадает под теорему 1. Тем не менее, и для него теорема оказалась верна. Полученные погрешности были аналогичны рис. 3. В частности, скорость сходимости для рассмотренных $m$ оказалась такой же. Количественная точность была несколько хуже рис. 3. Это было наиболее заметно для $m = 1.$ Для других $m$ ухудшение точности оказалось несущественным.
ЗАКЛЮЧЕНИЕ
Предложена специальная замена переменной интегрирования, которая кардинально повышает точность квадратуры средних. Для бесконечно гладких подынтегральных функций сходимость становится сверхстепенной. Для функций ограниченной гладкости закон сходимости остается степенным, но реализуется максимально достижимый порядок точности.
Проведем качественное сравнение предлагаемого подхода с другими методами повышения точности квадратур, перечисленными во введении.
Ни один из них не обеспечивает сверхстепенную сходимость. Поэтому для бесконечно гладких функций предлагаемый подход обеспечивает заведомо более высокую точность.
Использование поправок Эйлера–Маклорена требует большого количества априорной информации о подынтегральной функции. Нужно точно вычислять высокие производные и априори задавать число учитываемых поправок. Поэтому максимальный порядок точности реализуется, если известен класс гладкости подынтегральной функции и положение особенностей.
Напротив, предлагаемый подход единообразно применим к подынтегральным функциям как бесконечно гладким, так и имеющим ограниченную гладкость. Класс гладкости и положение особенностей заранее знать не нужно.
Повышать порядок точности с помощью ричардсоновской экстраполяции можно только на достаточно подробных сетках, на которых уже реализуется теоретическая сходимость, но еще не достигнуты ошибки округления. На грубых сетках применение экстраполяции может даже ухудшить точность.
В предлагаемом методе даже на грубых сетках наблюдается сходимость, причем достаточно быстрая. Количественное сравнение экстраполяции по Ричардсону и предлагаемого метода для функций ограниченной гладкости выходит за рамки данной работы.
Работа поддержана грантом Программой стратегического академического лидерства РУДН.
Список литературы
Калиткин Н.Н., Альшина Е.А. Численные методы. Т. 1. Численный анализ. М.: Академия, 2013.
Калиткин Н.Н., Альшин А.Б., Альшина Е.А., Рогов Б.В. Вычисления на квазиравномерных сетках. М.: Физматлит, 2005.
Trefethen L.N., Weideman J.A.C. // SIAM Rev. E. 2014. V. 56. No. 3. P. 385.
Kalitkin N.N., Kolganov S.A. // Dokl. Math. 2017. V. 95. No. 2. P. 157.
Kalitkin N.N., Kolganov S.A. // Math. Models Comp. Simul. 2018. V. 10. No. 4. P. 472.
Sag T.W., Szekeres G. // Math. Comp. 1964. V. 18. P. 245.
Sidi A. // Intern. Ser. Numer. Math. V. 112. P. 359.
Iri M., Moriguti S., Takasawa Y. // J. Comp. Appl. Math. 1987. V. 17. P. 3.
Mori M. // Publ. Res. Inst. Math. Sci. Kyoto Univ. 1978. V. 14. P. 713.
Белов А.А., Калиткин Н.Н., Хохлачев В.С. // Препринты ИПМ им. М.В. Келдыша. 2020. № 75.
Хохлачев В.С., Белов А.А., Калиткин Н.Н. // Изв. РАН. Сер. физ. 2021. Т. 85. № 2. С. 282.
Дополнительные материалы отсутствуют.
Инструменты
Известия РАН. Серия физическая
Пн.-пт.: 10.00-19.00