Повышение точности экспоненциально сходящихся квадратур

Аннотация

Вычисление одномерных интегралов возникает во многих задачах физики и техники. Для этого чаще всего используются простейшие квадратуры средних, трапеций и Симпсона на равномерной сетке. Для интегралов от периодических функций по полному периоду сходимость этих квадратур резко ускоряется и зависит от числа шагов сетки по экспоненциальному закону. В данной работе получены новые асимптотически точные оценки погрешности таких квадратур. Они учитывают расположение и кратность полюсов подынтегральной функции в комплексной плоскости. Построено обобщение этих оценок на случай, когда априорная информация о полюсах подынтегральной функции отсутствует. Описана процедура экстраполяции погрешности, которая кардинально ускоряет сходимость квадратур. Библ. 19. Фиг. 3.

Толық мәтін

1. Введение

1.1. Прикладные задачи

Рассмотрим задачу о вычислении однократного интеграла от функции f (x) по вещественному отрезку $x\in \left[\mathrm{0,}T\right]$

$I={\int }_0^{T}f\left(x\right)dx.$ (1)

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

1. Вычисление специальных функций математической физики: функции Ферми–Дирака, равные моментам фермиевского распределения, гамма-функция, цилиндрические функции, — W-функции Ламберта и ряд других.
2. Расчет скорости реакции по заданному сечению.
3. Расчет коэффициентов Фурье, преобразований Фурье и Лапласа заданной функции.
4. Численное решение интегральных уравнений, как корректно поставленных, так и некорректных.
5. Решение краевых задач для уравнений в частных производных (включая задачи на собственные значения), записанных в форме интегральных законов сохранения и т.д.

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

1.2. Вычисление квадратур

Как правило, для сеточного вычисления квадратур используют методы трапеций, средних и Симпсона на равномерной сетке [3]. Для этих методов хорошо известна мажорантная оценка погрешности. Для непериодических C²-гладких функций погрешность формул трапеций и средних есть O(N^–2), для формулы Симпсона — O(N^–4), если подынтегральная функция является C ⁴-гладкой. Здесь N — число шагов сетки. Если f (x) можно периодически продолжить за пределы отрезка (0,T ), и это периодическое продолжение непрерывно вместе с p производными, то порядок точности указанных квадратур повышается на p единиц.

Если подынтегральная функция является периодической (т.е. периодическое продолжение является бесконечно гладким), и интеграл вычисляется по полному периоду, то сходимость резко ускоряется. Зависимость погрешности от шага становится не степенной, а экспоненциальной $~exp(-const · N)$ [5, 6, 18]. Это означает, что при уменьшении шага вдвое число верных знаков в ответе примерно удваивается. Эта скорость сходимость аналогична таковой в методе Ньютона.

Аналогичной скоростью сходимости обладают спектральные методы, основанные на аппроксимации первообразной полиномами Чебышёва I рода (см., например, [7–13, 15–17]). В методах этого класса для первообразной составляют дифференциальное уравнение, решение которого представляют в виде разложения по базису указанных полиномов с некоторыми коэффициентами. Однако трудоемкость расчета оказывается довольно высокой. Кроме того, система уравнений относительно коэффициентов порой оказывается плохо обусловленной, поэтому для ее обращения используют регуляризацию.

1.3. Контроль точности и экстраполяция погрешности

1. Проведем расчеты на нескольких сетках с числом шагов N, 2N, 4N…, т.е. с шагом h, h/2, h/4… Получим последовательность сеточных значений квадратур I_N, I_2N, I_4N… Если, начиная с некоторого достаточно большого N, эта последовательность стремится к некоторому пределу, то можно констатировать сходимость квадратуры.

Если известна зависимость погрешности $\delta$ от числа шагов, то можно построить асимптотически точные оценки точности. Например, для разностных формул порядка точности p имеем $\delta ~ N^{-p}$ т.е. при уменьшении шага вдвое главная асимптотика погрешности убывает в 2^p раз. Отсюда следуют хорошо известные оценки погрешности по методу Ричардсона [4, 14]

${\delta }_{2N}=\frac{I_{2N}-I_N}{2^p-1}$. (2)

либо

${\delta }_N=2^p\frac{I_{2N}-I_N}{2^p-1}.$ (3)

Форма записи (3) более предпочтительна, поскольку она работоспособна не только на участке регулярной сходимости $\delta ~ N^{-p}$, но и по достижении ошибок округления.

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

${\tilde{I}}_{2N}=I_{2N}+{\delta }_{2N}.$ (4)

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

2. Если известны производные подынтегральной функции в точках x = 0 и x = T порядка до 2M -1 включительно, то погрешность можно оценить по формуле Эйлера-Маклорена

$\delta =\sum _{m=1}^M{c}_m{h}^{2m}\left(f^{(2m-1)}\left(T\right)-f^{(2m-1)}\left(0\right)\right).$ (5)

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

3. Для случая, когда квадратуры средних и трапеций сходятся экспоненциально, известны оценки точности Трефетена–Вайдемана [18]. Однако они не являются асимптотически точными и поэтому непригодны для экстраполяции погрешности.

Нами ранее были построены апостериорные асимптотически точные оценки для сеточных методов с экспоненциальной сходимостью [19]. Они напоминают таковые в методе Ричардсона. Эти оценки имеют вид

${\delta }_N=I_{2N}-I_N,$ (6)

${\delta }_{4N}=\frac{{(I_{4N}-I_{2N})}^3}{{(I_{2N}-I_N)}^2}.$ (7)

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

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

2. Погрешность квадратуры трапеций

Пусть подынтегральная функция f (x) является периодической на отрезке (0,T). С помощью замены переменной $z = exp(2\pi ix / T)$ перейдем от (1) к интегралу по единичной окружности на комплексной плоскости z

$I={\int }_{\left|z\right|=1}g\left(z\right)dz,$ (8)

где

$g\left(z\right)=\frac{T}{2\pi iz}f\left(\ln z\right).$(9)

Введем на окружности $\left|z\right|=1$ сетку $z_n=\exp \left(2\pi in/N\right),$ $n=\mathrm{0,1,}\dots ,N-\mathrm{1,}$  шаг сетки равен $h_n=2\pi i{N}^{-1}{z}_n.$ Запишем квадратуру трапеций для интеграла (8)

$I_N=2\pi i{N}^{-1}{\sum }_{n=0}^{n=N-1}g\left(z_n\right)z_n,$ (10)

2.1. Оценка Трефетена–Вайдемана

В работе [18] доказано следующее утверждение. Пусть g(z) имеет полюс первого порядка при z = 0, функция zg(z) аналитична и удовлетворяет неравенству |zg (z )| < M, M = const в кольце R^-1 < |z| < R для некоторого R. Тогда для квадратуры (10) справедлива мажорантная оценка погрешности

${\delta }_N=\left|I-I_N\right|\le 2M{\left(R^N-1\right)}^{-1}.$ (11) 

Наш анализ показал следующее [2]. Во-первых, зависимость d_N от N асимптотически точна, если все конечные особые точки g(z) суть полюсы первого порядка, и ближайшая к контуру особая точка попадает на границу указанного кольца. В этом случае константа M формально обращается в бесконечность. Во-вторых, если ближайшая к контуру особая точка является полюсом порядка 2 и выше, то оценка (11) теряет применимость.

2.2. Единственный полюс

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

Сначала рассмотрим простейший случай: пусть g(z) имеет единственный полюс порядка q, расположенный в точке z = a внутри единичной окружности, и не имеет других конечных особых точек. Отметим, что при q ³ 2 даже этот простейший случай выходит за пределы применимости оценки (11).

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

$g\left(z\right)=f\left(z\right){\left(z-a\right)}^{-q}={\left(\left(q-1\right)!\right)}^{-1}{\partial }_a^{\left(q-1\right)}\left(f\left(z\right){\left(z-a\right)}^{-1}\right),$(12)

где функция f(z) не имеет конечных особых точек. Тогда нетрудно найти точное значение интеграла $I={\left(\left(q-1\right)!\right)}^{-1}{f}^{\left(q-1\right)}\left(a\right).$

Найдем численное значение интеграла. Разложим в ряды по степеням z функции f(z) и (z – a)^–1

$f\left(z\right)={\sum }_{q=0}^{+\infty }{f}_q{z}^q,{\left(z-a\right)}^{-1}={\sum }_{k=0}^{+\infty }{a}^k{z}^{-k-1}.$(13)

Здесь  — коэффициенты ряда Тейлора функции f(z). Подставим в (10) представление (12), (13)

$I_N=2\pi i{N}^{-1}{\left(\left(q-1\right)!\right)}^{-1}{\partial }_a^{\left(q-1\right)}{\sum }_{q,k=0}^{+\infty }{f}_q{a}^k{\sum }_{n=0}^{N-1}{z}_n^{q-k}.$ (14)

Нетрудно показать, что сумма по n равна N, если g - k кратно N, и нулю в остальных случаях. Поэтому в двойной сумме по q,k остаются только те слагаемые, у которых q − k = sN где s — целое число. Его удобно выбрать в качестве нового индекса суммирования. Диапазон его изменения есть $-\infty < s < [q / N ],$, где квадратные скобки обозначают взятие целой части. Имеем

$I_N={\left(\left(q-1\right)!\right)}^{-1}{\partial }_a^{\left(q-1\right)}{\sum }_{q=0}^{+\infty }{f}_q{a}^q{\sum }_{s=-\infty }^{\left[q/N\right]}{a}^{-sN}.$ (15)

Внутренний ряд есть сумма геометрической прогрессии. Внешняя сумма в главной по N асимптотике равна ряду Тейлора для f (a). Поэтому главная по N асимптотика для  равна

${\left.I_N\right|}_{N\to \infty }={\left(\left(q-1\right)!\right)}^{-1}{\partial }_a^{\left(q-1\right)}\left(f\left(a\right){\left(1-a^N\right)}^{-1}\right).$ (16)

Найдем главную асимптотику погрешности квадратуры (10)

${\delta }_N={\left.I-I_N\right|}_{N\to \infty }={\left(\left(q-1\right)!\right)}^{-1}{\partial }_a^{\left(q-1\right)}\left(f\left(a\right)\left(1-{\left(1-a^N\right)}^{-1}\right)\right).$ (17)

2.3. Множественные полюсы произвольной кратности

Формула (17) обобщена на случай, когда подынтегральная функция имеет произвольное конечное число полюсов целого порядка как внутри, так и вне контура интегрирования. Пусть функция g(z) имеет полюсы в точках $a_{1,} a_{2,...,} a_m$ порядков $q_{1,} q_{2,...,} q_m$ соответственно. Тогда

$g\left(z\right)=f\left(z\right){\left(z-a_1\right)}^{-q_1}{\left(z-a_2\right)}^{-q_2}\dots {\left(z-a_m\right)}^{-q_m}.$ (18)

Здесь функция f (z) не имеет конечных особых точек.

Опуская громоздкие выкладки, приведем окончательный результат. Введем обозначения $q = q_1 + q_2 +\dots + q_m , A_{ j}= (a_j -a_1)(a_j -a_2 )\dots (a_j -a_j-1)(a_j -a_{j+1})\dots (a_j -a_m )$ и

$L={\prod }_{s=1}^m\frac{1}{\left(q_s-1\right)!}\frac{{\partial }^{q-m}}{\partial a_1^{q_1-1}\partial a_2^{q_2-1}\dots \partial a_m^{q_m-1}}.$(19)

Главная по  асимптотика погрешности квадратуры (10) от функции (18) имеет вид

${\left.{\delta }_N\right|}_{N\to \infty }=L{\sum }_{k}f\left(a_k\right)A_k^{-1}\left(1-{\left(1-a_k^N\right)}^{-1}\right)-L{\sum }_{j}f\left(a_j\right)A_j^{-1}{\left(a_j^N-1\right)}^{-1}.$ (20)

Здесь первая сумма берется по всем индексам k таким, что |a_k| < 1; вторая сумма — по индексам j, для которых |a_j| >1.

По построению, оценка (20) не мажорантна, а асимптотически точна. Она применима к существенно более широкому классу подынтегральных функций, чем оценка (11) из работы [18]. Поэтому оценка (20) кардинально превосходит мировой уровень.

2.4. Неизвестные параметры

Чтобы применять оценки (17) и (20), нужно знать положения a_k и кратности q_k полюсов подынтегральной функции. Такая априорная информация известна не всегда. Например, если подынтегральная функция сама является решением некоторой сеточной задачи, то отыскание ее полюсов в комплексной плоскости представляет самостоятельную проблему (причем более сложную, чем расчет исходного интеграла).

Однако зависимость d от N в оценке (17) является асимптотически точной. Это позволяет рассчитать неизвестное a_k для полюса, наиболее близкого к контуру интегрирования.

Пусть проведен расчет на наборе сеток с числами узлов N, 2N, 4N,… Обозначим через I_end значение квадратуры с самой последней сетки. В силу исключительно быстрой сходимости это значение намного точнее расчетов с предыдущих сеток. Поэтому точность квадратур, вычисленных на всех сетках, кроме последней, можно оценить как

${\tilde{\delta }}_N=I_{end}-I_N.$ (21)

Построим оценку точности для I_end.

Простой полюс. Пусть подынтегральная функция в (8) имеет один простой полюс в точке a внутри единичной окружности. Случай полюса, расположенного вне единичной окружности |a| > 1 сводится к предыдущему заменой $a \to 1 / a.$  Поэтому эти случаи можно рассматривать вместе.

Согласно (17), для погрешности квадратуры (10) имеет место оценка

${\left.{\delta }_N\right|}_{N\to \infty }={\left.\left(I-I_N\right)\right|}_{N\to \infty }=f\left(a\right)\left(1-{\left(1-a^N\right)}^{-1}\right).$ (22)

С другой стороны, величину d_N можно апостериорно оценить согласно (21).

Удобно ввести обозначения A = f (a), t = a^N. Обе эти величины неизвестны. Запишем равенства (21), (22) для двух последовательных сеток, содержащих N и 2N шагов

$$\begin{gathered} {\tilde{\delta }}_N=A\left(1-{\left(1-t\right)}^{-1}\right), \\ {\tilde{\delta }}_{2N}=A\left(1-{\left(1-t^2\right)}^{-1}\right). \end{gathered}$$ (23)

Равенства (23) есть система уравнений относительно A, t. Ее решение имеет вид

$t=1/\left(\xi -1\right),A={\tilde{\delta }}_N\left(t-1\right)/t.$(24)

Здесь для удобства введено обозначение $\xi ={\tilde{\delta }}_N/{\tilde{\delta }}_{2N}$

Пользуясь найденными A и t, вычислим оценку погрешности, соответствующую сетке с 4N шагами

${\delta }_{4N}=A\left(1-{\left(1-t^4\right)}^{-1}\right)=A{t}^4/\left(t^4-1\right).$(25)

Формулу (25) можно использовать для контроля точности расчета на всех сетках, включая последнюю.

Полюс второго порядка. Аналогично рассмотрим случай полюса второго порядка. Пусть нам известны оценки погрешности на трех соседних сетках (содержащих N, 2N, 4N шагов), и нужно оценить погрешность на следующей сетке (8N шагов). Оценка (20) в этом случае примет вид

${\delta }_N=f^{\text{'}}\left(a\right)\left(1-{\left(1-a^N\right)}^{-1}\right)-f\left(a\right)N{a}^{N-1}/{\left(1-a^N\right)}^2.$ (26)

С другой стороны, для d_N справедлива апостериорная оценка (21).

Введем обозначения $A = f \text{'}\left(a\right), B = f \left(a\right) / a, t = 1 / (1-a^N ).$ Эти величины являются неизвестными. Запишем оценки (21), (26) для сеток N, 2N, 4N

$$\begin{gathered} -{\delta }_N/\left(t-1\right)=A+BNt, \\ -{\delta }_{2N}\left(t^2-{\left(t-1\right)}^2\right)/{\left(t-1\right)}^2=A+2BN{t}^2/\left(t^2-{\left(t-1\right)}^2\right), \\ -{\delta }_{4N}\left(t^4-{\left(t-1\right)}^4\right)/{\left(t-1\right)}^4=A+4BN{t}^4/\left(t^4-{\left(t-1\right)}^4\right). \end{gathered}$$(27)

Равенства (27) есть система уравнений для определения A, B, t. Исключим из нее A и B и после несложных, но громоздких выкладок получим уравнение относительно t

$E\left(t\right)\left(t-1\right)\left(\xi t^2+\left(\xi -\eta \right)t+\xi \right)+{\delta }_{2N}\left(t-1\right)\left(t+1\right)+{\delta }_{4N}=\mathrm{0,}$(28)

где $\xi =4{\delta }_{4N}+2{\delta }_N-6{\delta }_{2N}$ и $\eta =8{\delta }_{4N}+6{\delta }_N-16{\delta }_{2N}.$

Уравнение (28) допускает точное решение в радикалах. Однако корни очень громоздки, поэтому мы их не приводим. Для их вычисления мы разработали программу на языке Matlab/Octave [1].

Альтернативным способом решения уравнения (28) является метод Ньютона. В нем важным аспектом является выбор начального приближения. Из явного вида уравнения (28) следует, что все 4 корня “слипаются” в один кратный корень t = 1, при $N\to \infty$ в силу сходимости метода средних. Поэтому в качестве начального приближения целесообразно брать $t=1.$ Однако численное отыскание нужного корня уравнения $E\left(t\right)=0$ сталкивается с рядом трудностей, решение которых выходит за рамки данной работы.

Наконец, по найденным значениям A, B, t, вычислим оценку точности по формуле (26), заменяя $N\to 8N.$ Эту оценку точности можно использовать для всех сеток, включая последнюю.

Полюсы высокой кратности. Аналогично, но существенно более громоздко рассматриваются случаи, когда ближайшие полюсы имеют кратность 3 и выше. Например, для полюса 3-го порядка оценка погрешности содержит члены $~a^{-N},$ $N{a}^{-N},$ $N^2{a}^{-N}.$ Система уравнений, аналогичная (23), (27), включает 4 уравнения. После исключения части неизвестных она приводит к алгебраическому уравнению 6-й степени. Для его решения целесообразно использовать метод Ньютона. Соответствующие формулы весьма громоздки, поэтому мы их не приводим.

3. Формула средних

Построенные выше оценки непосредственно переносятся на формулу средних. В самом деле, подынтегральная функция является периодической на единичной окружности. Поэтому поворот комплексной плоскости на произвольный угол не изменяет значение интеграла. В этом случае формула трапеций переходит в формулу средних при повороте на угол $\pi /N.$ Это эквивалентно “переносу” полюсов подынтегральной функции $a_q\to a_q{e}^{-i\pi /N}.$

4. Экстраполяция точности

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

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

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

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

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

Отметим один нюанс. Алгоритм, описанный в п. 2.4, есть некоторый сеточный метод расчета a^N. Входными данными этого алгоритма являются апостериорные оценки ${\tilde{\delta }}_N,$ которые вычисляются с ошибкой. Поэтому найденные a^N являются не точными, а приближенными. Поэтому экстраполированная квадратура оказывается менее точной, чем при использовании оценки (20) с точно известными положениями полюсов.

5. Примеры

Рассмотрим вычисление тестового интеграла

$I={\left(2\pi i\right)}^{-1}{\int }_{\left|z\right|=1}th\left({\left(z-b\right)}^{-1}\right){\left(z-a\right)}^{-q}dz,$ (30)

при различных значениях параметров a, b и q.

Простой полюс. Пусть подынтегральная функция в (30) имеет один простой полюс (q = 1) в точке $a={\pi }^{-1}\approx 0.3183\dots$ и существенно особую точку (СОТ) $b={10}^3{\pi }^{-1}\approx 318.3098\dots .$ Таким образом, ближайшей к контуру интегрирования $\left(\left|z\right|=1\right)$ особой точкой является указанный полюс. Интеграл (30) в данном случае имеет точный ответ. Он равен $I=th\left({\left(a-b\right)}^{-1}\right)\approx -0.003144727024527\dots .$ На фиг. 1 показана погрешность вычисления этого интеграла в зависимости от числа N шагов сетки. Масштаб графика полулогарифмический; экспоненциальной сходимости соответствует прямая линия.

Фиг. 1. Погрешность расчета тестового интеграла (30) при $a={\pi }^{-1},$ $b={10}^{-3}{\pi }^{-1}$ и $q=1.$  Обозначения — см. текст.

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

Темными кружками показана оценка (11) для $R=2.8$ и, соответственно, $M\approx 0.00316946\dots .$ Эта линия является прямой. Видно, что с ростом N оценка убывает медленнее фактической погрешности. При этом оценка (11) мажорирует фактическую погрешность.

Светлыми квадратами показана погрешность квадратуры (10) после апостериорной экстраполяции (29). Видно, что по точности она кардинально превосходит исходную квадратуру трапеций. В пределах графика выигрыш составляет от $~{10}^3$ до $~{10}^{18}$ раз. При дальнейшем увеличении N преимущество еще увеличивается. Погрешность порядка ошибок округления двойной точности $~{10}^{-16}$ достигается уже при $N~18.$

Темными квадратами показана погрешность квадратуры (10) после априорной экстраполяции (29) с помощью оценки (20). Она тоже представляется прямой линией. Такая квадратура сходится кардинально быстрее остальных. Точность $~{10}^{-16}$ достигается уже при $N~7.$

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

Полюс второго порядка. Рассмотрим тестовый интеграл (30) с полюсом порядка q = 2 в точке $a={\pi }^{-1}$ и СОТ в точке $b={10}^3{\pi }^{-1}.$ Его точное значение равно $I=-{\cosh }^{-2}\left({\left(a-b\right)}^{-1}\right){\left(a-b\right)}^{-2}\approx -0.000009889275459\dots .$ Соответствующие погрешности показаны на фиг. 2. Обозначения соответствуют фиг. 1. Масштаб графика также полулогарифмический.

Фиг. 2. Погрешности расчета тестового интеграла (30) при $a={\pi }^{-1},$ $b={10}^{-3}{\pi }^{-1}$ и $q=2.$  Обозначения — см. текст.

Из фиг. 2 видно, что поведение всех зависимостей на графике, кроме оценки (11), аналогично предыдущему случаю (фиг. 1). Кривая, соответствующая оценке (11), пересекает кривую погрешности относительно точного значения интеграла. Это означает, что оценка (11) в данном случае не является мажорантной.

Видно, что процедура апостериорной экстраполяции дает значительный выигрыш в точности: от $~{10}^7$ до $~{10}^{20}$ раз. Таким образом, экстраполяционная квадратура (29) кардинально превосходит обычную квадратуру трапеций (10) по точности. Видно также, что априорная экстраполяция дает кардинальное ускорение сходимости. Выигрыш по точности составляет десятки порядков.

Два простых полюса. Для иллюстрации работы алгоритма рекуррентной экстраполяции рассмотрим следующий пример

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

Положим $a_1={\pi }^{-1}\approx 0.3183\dots ,$ $a_2=10\pi \approx 31.4159\dots$ и $b={10}^3{\pi }^{-1}\approx 318.3098\dots .$ Точный ответ в этом случае будет равен $I=th\left({\left(a_1-b\right)}^{-1}\right){\left(a_1-a_2\right)}^{-1}\approx 0.000101124374255\dots .$

Погрешности, полученные в этом примере, представлены на фиг. 3. Масштаб графика полулогарифмический. Видно, что оценка (11) является мажорантной, что согласуется с теорией. Апостериорная экстраполяция проводилась реккурентно: первым исключался ближайший полюс, затем второй. Из фиг. 3 видно, что первая экстраполяция погрешности (светлые квадраты) дает существенное уточнение квадратуры по сравнению с фактической погрешностью. Вторая экстраполяция (пунктирная линия) также уточняет расчет, но заметно слабее. Как и в предыдущих примерах, хорошо видно, что априорная экстраполяция дает кардинальный выигрыш по точности (десятки порядков) и значительно превосходит апостериорную.

Фиг. 3. Погрешности расчета тестового интеграла (31) при $a_1={\pi }^{-1},$ $a_2=10\pi ,$ $b={10}^3{\pi }^{-1}.$ Обозначения — см. текст.

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

6. Заключение

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

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

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