Журнал вычислительной математики и математической физики, 2020, T. 60, № 4, стр. 553-566

1. ВВЕДЕНИЕ

С появлением фундаментальных исследований К.И. Бабенко [1], [2] проблематика численного решения эллиптических краевых задач вступила в ту стадию своего естественного развития, когда пренебрежение ее высокой теоретической оснащенностью (наличием шаудеровских оценок гладкости) становится тормозом дальнейшего развития, и оценка практической эффективности численных алгоритмов уже невозможна без использования для этой цели понятий александровского $m$-поперечника ${{\alpha }_{m}}(X)$ и колмогоровской $\varepsilon $-энтропии ${{H}_{\varepsilon }}(X)$ компакта $X$ решений этих задач.

Наилучшее финитное описание функционального компакта $X$ ассоциируется всегда с асимптотиками двух числовых параметров ${{\alpha }_{m}}(X)$ и ${{H}_{\varepsilon }}(X)$, характеризующими соответственно степень $m$-мерности $X$ и минимальный в битах объем информации, требуемый для различения элементов в $X$ с точностью $\varepsilon > 0$. При этом убывание ${{\alpha }_{m}}(X)$ к нулю при $m \to \infty $ тем быстрее, а рост ${{H}_{\varepsilon }}(X)$ к бесконечности при $\varepsilon \to 0$ тем медленнее, чем выше запас гладкости $X$; для компакта $X$ бесконечной гладкости ${{\alpha }_{m}}(X)$ убывает к нулю экспоненциально, а рост ${{H}_{\varepsilon }}(X)$ к бесконечности определяется фиксированной степенью $log(1{\text{/}}\varepsilon )$ [2], [3]. При этом качество финитизаций $X$ характеризуется сопоставлением точности $\varepsilon > 0$ и объема числовых данных с асимптотиками параметров ${{\alpha }_{m}}(X)$ и ${{H}_{\varepsilon }}(X)$ – аппроксимационными возможностями компакта $X$.

Смысловые потенции понятий ${{\alpha }_{m}}(X)$ и ${{H}_{\varepsilon }}(X)$, сыграв решающую роль в становлении точных представлений о предельных возможностях численных алгоритмов, открывают доступ к принципиально новому – ненасыщаемому – типу вычислительных алгоритмов, до недавнего времени не входившему в арсенал вычислительной практики. Сообщение К.И. Бабенко об этом открытии появилось в 1975 г. (см. Докл. АН СССР, 1975, Т. 221. № 1).

Краткое изложение основ теории ненасыщаемых численных методов (и алгоритмов) содержится в книге [2], второе издание которой появилось в 2002 г. Отличительная черта ненасыщаемых численных методов – отсутствие главного члена погрешности, и как результат – способность автоматически подстраиваться к любым аппроксимационным возможностям компакта $X$ решений задач, определяемым асимптотикой убывания ${{\alpha }_{m}}(X)$ к нулю (поперечник ${{\alpha }_{m}}(X)$ определяется как нижняя грань $\varepsilon $-сдвигов компакта $X$ в компакт размерности, не большей $m$ [1], [2]). Скорость убывания ${{\alpha }_{m}}(X)$ к нулю сравнивается при этом с числом $m$ свободных параметров конечномерного описания $X$: она тем выше, чем бóльший “запас” гладкости имеет $X$. Так что с ростом гладкости $X$ (при прочих равных условиях) ненасыщаемый численный метод самосовершенствуется, черпая приращение своей практической эффективности непосредственно в дифференциальной природе $X$ (феномен ненасыщаемости [4]). В результате избыточная (экстраординарная) гладкость $X$, прежде находившаяся на периферии насущных интересов вычислений, становится активным персонажем. Причем пик эффективности – экспоненциальная сходимость (или сверхсходимость) – достигается на компактах бесконечно гладких $X$. И это принципиально отличает ненасыщаемые численные методы от методов с главным членом погрешности, т.е. насыщаемых. В итоге, открывая новые перспективы для вычислительной практики, ненасыщаемые численные методы создают серьезную основу для построения компьютерных числовых ответов гарантированного качества [5].

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

К сказанному следует добавить, что, поскольку операция замены вещественного числа конечной дробью (округление) является вынужденной и необходимой, и без потери информации не обходится, переход к цифровым вычислениям оказывается моментом небезобидным. И потому в качестве составной части он обязан содержать компонент “хорошей обусловленности” (или корректности) компьютерной постановки задачи, определяемый устойчивостью к ошибкам округлений и характеризуемый параметром, называемым числом обусловленности [5]. Если это число велико, то от такого способа финитизации задачи следует отказаться, потому как конструируемый числовой ответ не может внушать доверия из-за плохой устойчивости алгоритма.

В работе на примере численного решения задачи Неймана указана методика построения ненасыщаемых алгоритмов численного решения эллиптических краевых задач для уравнения Лапласа в ${{C}^{\infty }}$-гладких осесимметричных областях достаточно произвольной формы. В общих чертах методика описана в работах [6], [7]. В ее основе – методы теории гармонического потенциала [8], позволяющие сводить задачи к эквивалентным им интегральным уравнениям, с последующим отысканием их решений. Предлагаемый подход в отличие от конечно-разностного метода, связанного с аппроксимацией дифференциального (неограниченного) оператора, имеет дело с аппроксимацией оператора интегрального, т.е. ограниченного. А это в сумме с ненасыщаемой дискретизацией последнего приводит к хорошо обусловленной системе линейных алгебраических уравнений, которая затем эффективно решается быстросходящимся итерационным методом. При этом не возникает проблем с решением систем “большой” размерности, а стало быть, и коллизий точности конструируемого числового ответа и объема перерабатываемой битовой информации, если решение задачи достаточно гладкое, например, бесконечно дифференцируемое [3].

Имеющиеся в настоящее время случаи применения ненасыщаемых численных методов в конкретных задачах [6], [9], [10], хотя и немногочисленны, но убеждают в их содержательности и высокой эффективности на практике.

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

2. ПОСТАНОВКА ПРОБЛЕМЫ

Пусть $x,\xi \in {{\mathbb{R}}^{3}}$, $x = (x,y,z)$, $\xi = (\xi ,\eta ,\zeta )$; $\omega \subset {{\mathbb{R}}^{3}}$ – связная область с осью симметрии $z$, ограниченная ${{C}^{\infty }}$-гладкой замкнутой поверхностью вращения $\partial \omega $; величины $r = \sqrt {{{x}^{2}} + {{y}^{2}}} $ и $z$ – инварианты группы вращения $\omega $ относительно оси $z$. Меридиональное сечение $\partial \omega $ – параметризованная кривая $\gamma :[0,\pi ] \to {\text{\{ }}r(s),z(s){\text{\} }}$, $r \geqslant 0$, $dz{\text{/}}ds \geqslant 0$ и $\gamma (s) \in {{C}^{\infty }}[0,\pi ]$. Точки $\gamma (0)$ и $\gamma (\pi )$ суть полюсы $\partial \omega $, и функции $r(s)$, $z(s)$ имеют $2\pi $-периодические ${{C}^{\infty }}$-гладкие продолжения (нечетное и четное соответственно) с отрезка $[0,\pi ]$ на $[0,2\pi ]$. Положение точек $x = (r,z)$ и $\xi = (\rho ,\zeta )$ на $\gamma $ определим координатами $s$ и $\sigma $: $r = r(s)$, $z = z(s)$; $\rho = r(\sigma )$, $\zeta = z(\sigma )$ и $0 \leqslant s$, $\sigma \leqslant \pi $. Введем также обозначения: $\rho {\text{'}} \equiv d\rho {\text{/}}d\sigma $, $\zeta {\text{'}} \equiv d\zeta {\text{/}}d\sigma $, $\delta = \sqrt {\rho {{{\text{'}}}^{2}} + \mathop {\zeta {\kern 1pt} }\nolimits^ {{{\text{'}}}^{2}}} $; $r{\text{'}} \equiv dr{\text{/}}ds$, $z{\text{'}} \equiv dz{\text{/}}ds$, $d = \sqrt {r{\kern 1pt} {{{\text{'}}}^{2}} + z{\kern 1pt} {{{\text{'}}}^{2}}} $. Нормаль $N$ на $\partial \omega $ считается направленной в ${{\mathbb{R}}^{3}}{\backslash }\omega $.

Пусть $F(x)$ – непрерывная функция вне области $\omega $. Прямое значение $F(x)$ на $\partial \omega $ (если оно существует) обозначим через $\bar {F}(x) = {{\left. {F(x)} \right|}_{{x \in \partial \omega }}}$. Для предельного извне $\partial \omega $ и прямого на $\partial \omega $ значений производной $\partial F{\text{/}}\partial N = N \times {{\nabla }_{x}}F$ (если они существуют) примем обозначения: $({{N}_{ - }}\bar {F})(x)$ и $(N\bar {F})(x)$.

Решение внешней задачи Неймана для уравнения Лапласа в ${{C}^{\infty }}$-гладкой осесимметричной области $\omega $ с данными ${{\left. {{{N}_{ - }}\bar {\varphi }} \right|}_{{x \in \partial \omega }}} = f(x)$ и условием $\varphi \to 0$ при $\left| x \right| \to \infty $ будем отыскивать в виде потенциала простого слоя $\varphi (x) = V[\chi ](x)$ с достаточно гладкой плотностью $\chi (\xi )$, инвариантной относительно группы вращения $\omega $. Решение задачи существует и единственно. Гладкость решения $\varphi (x)$ определяется гладкостью функции $f(x)$. Обычно полагают [8], что $f(x) \in {{C}^{\alpha }}(\partial \omega )$, т.е. – пространству непрерывных функций, удовлетворяющих на $\partial \omega $ условию Гёльдера с показателем $0 < \alpha < 1$.

Пусть $f(s) \equiv f(r(s),z(s))$ – достаточно гладкая $2\pi $-периодическая четная функция, тогда интегральное уравнение для внешней задачи Неймана имеет вид [8]

Здесь оператор $N\bar {V} \equiv K$ – интегральный, который согласно инвариантности относительно группы вращения области $\omega $, задается равенством Функции $E(q)$ и $D(q)$ – полные эллиптические интегралы с модулем $q \in [0,1]$, где Эффективные алгоритмы вычисления функций $E(q)$ и $D(q)$ указаны в работе [11].

Приняты также следующие обозначения:

3. ДИСКРЕТИЗАЦИЯ УРАВНЕНИЯ (2.1) И ФОРМУЛИРОВКА РЕЗУЛЬТАТА

Пусть $C \equiv C[0,2\pi ]$ – класс вещественных $2\pi $-периодических непрерывных функций с чебышëвской нормой $\left\| {\, \cdot \,} \right\|$, а ${{C}_{ + }} \equiv {{C}_{ + }}[0,\pi ] \subset C$ – множество четных функций. Пусть ${{\mathcal{T}}^{m}} \subset {{C}_{ + }}$ – подпространство тригонометрических многочленов, порядок которых не выше $m$, $m \geqslant 0$ – целое, ${{Q}_{m}}:{{C}_{ + }} \to {{\mathcal{T}}^{m}}$ – проектор, ${{e}_{m}}(g) = \mathop {inf}\limits_{{{T}_{m}} \in {{\mathcal{T}}^{m}}} \left\| {g - {{T}_{m}}} \right\|$ – наилучшее (чебышёвское) приближение функции $g$ из ${{C}_{ + }}$ многочленами.

Ядро $k(s,\sigma )$ – слабосингулярно и потому оператор $K:{{C}_{ + }} \to {{C}_{ + }}$ компактен и его норма вычисляется по формуле

Дискретизацию уравнения (2.1) осуществим по следующей схеме. Зададим узлы

и линейное отображение Рассмотрим интерполяционный многочлен Лагранжа функции $g$ из ${{C}_{ + }}$: Здесь – есть ядро Дирихле; ${{w}_{k}}({{s}_{j}}) = {{\delta }_{{kj}}}$, $0 \leqslant k$, $j \leqslant m$.

Многочлен ${{Q}_{m}}(s;Jg)$ определяет проектор ${{Q}_{m}}:{{C}_{ + }} \to {{\mathcal{T}}^{m}}$, а $\left\| {{{Q}_{m}}} \right\|$ – его норма. Ясно, что ${{Q}_{m}}{{w}_{k}} \equiv {{w}_{k}}$ и $J{{Q}_{m}}g \equiv Jg$.

Согласно неравенству Лебега [2] имеем оценку погрешности:

Далее, в силу линейности оператора $K:{{C}_{ + }} \to {{C}_{ + }}$, имеем

Пусть Применяя оператор $J$ к обеим частям уравнения (2.1), находим Здесь $A = ({{a}_{{ik}}})$ – матрица; ${{a}_{{ik}}} = K[{{w}_{k}}]({{s}_{i}})$ ($0 \leqslant k \leqslant m$, $0 \leqslant i \leqslant m$).

Матрица $A$ определяет линейный оператор $A:{{\mathbb{R}}^{{m + 1}}} \to {{\mathbb{R}}^{{m + 1}}}$, который рассматривается в качестве дискретизации $K:{{C}_{ + }} \to {{C}_{ + }}$; матрица $A$ полностью заполнена, в отличие от используемых в конечно-разностных методах сильно разреженных матриц.

Отбрасывая в (3.3) погрешность $\varrho \in {{\mathbb{R}}^{{m + 1}}}$ и, обозначая приближенное значение $J\chi \in {{\mathbb{R}}^{{m + 1}}}$ через $\bar {\chi } \in {{\mathbb{R}}^{{m + 1}}}$, получим дискретизацию интегрального уравнения (2.1):

Пусть ${{\left| B \right|}_{\infty }}$ – чебышёвская норма обратимой матрицы $B:{{\mathbb{R}}^{{m + 1}}} \to {{\mathbb{R}}^{{m + 1}}}$; меру ее обусловленности определим числом $\varkappa (B) = {{\left| B \right|}_{\infty }}{{\left| {{{B}^{{ - 1}}}} \right|}_{\infty }}$.

Уравнение (2.1) однозначно разрешимо, существует обратный оператор ${{(I + K)}^{{ - 1}}}:{{C}_{ + }} \to {{C}_{ + }}$, причем $\left\| {{{{(I + K)}}^{{ - 1}}}} \right\| \leqslant Q < \infty $ и выполняется оценка $\left\| \chi \right\| \leqslant Q\left\| f \right\|$.

В этом случае справедливы следующие нетривиальные результаты [6].

Теорема 1. Существуют постоянные ${{m}_{0}}$, ${{q}_{0}}$, ${{\varkappa }_{0}}$, ${{c}_{0}}$, такие что при $m \geqslant {{m}_{0}}$ существует матрица $\mathop {(I + A)}\nolimits^{ - 1} $ и имеют место следующие оценки

Постоянные ${{q}_{0}}$, ${{\varkappa }_{0}}$, ${{c}_{0}}$ не зависят от $m$ и эффективно вычисляются.

Следствие 1. Построенный метод численного решения задачи (2.1) ненасыщаем.

Доказательство. Считаем известными классические определения александровского ${{\alpha }_{m}}(X)$ и колмогоровского ${{\varkappa }_{m}}(X)$ $m$-поперечников функционального компакта $X$ и теорему о связи их величин при $m \to \infty $ между собой [12]: ${{c}_{1}}{{\varkappa }_{m}}(X) \leqslant {{\alpha }_{m}}(X) \leqslant {{c}_{2}}{{\varkappa }_{m}}(X)$ с не зависящими от $m$ постоянными ${{c}_{1}} > 0$, ${{c}_{2}} > 0$.

Упомянутая связь дает возможность использовать колмогоровский $m$-поперечник ${{\varkappa }_{m}}(X)$ в определении ненасыщаемости [2] вычислительного метода.

Действительно, для $k \geqslant 0$ и $0 < \alpha < 1$ рассмотрим класс четных $2\pi $-периодических $k$ раз непрерывно дифференцируемых функций $C_{ + }^{{k + \alpha }}(M)$, производные $k$-го порядка которых удовлетворяют условию Гёльдера с показателем $\alpha $ и константой $M$:

Тогда ${{X}^{k}} \equiv C_{ + }^{{k + \alpha }}(M) \cap {\text{\{ }}\psi \,|\,\psi \in {{C}_{ + }}[0,\pi ],\;\left\| \psi \right\| \leqslant 1{\text{\} }}$ – компакт в пространстве ${{C}_{ + }}[0,\pi ]$.

Пусть $\chi \in {{X}^{k}}$, ${{\varrho }_{m}}(\chi ) = \left\| {\chi (s) - {{Q}_{m}}(s;J\chi )} \right\|$ и ${{\epsilon }_{m}}({{X}^{k}}) = sup{\text{\{ }}{{\varrho }_{m}}(\chi )\,|\,\chi \in {{X}^{k}}{\text{\} }}$ – точность приближения элемента $\chi $ и точность метода на классе ${{X}^{k}}$ соответственно.

Компакты ${{X}^{k}}$, где $0 \leqslant k \leqslant \infty $, вложены друг в друга: если ${{k}_{1}} < {{k}_{2}}$, то ${{X}^{{{{k}_{2}}}}} \subset {{X}^{{{{k}_{1}}}}}$.

Ненасыщаемость метода решения уравнения (2.1) состоит в проверке двух условий:

Убедимся в этом. Согласно неравенству Джексона [2], из правой оценки в (3.5) следует С другой стороны, оценка ${{\varrho }_{m}}(\chi ) \geqslant {{e}_{m}}(\chi )$ и оценка снизу [12, c. 42] для ${{\varkappa }_{m}}({{X}^{k}})$ дают (постоянные ${{c}_{3}} > 0$, ${{b}_{k}} > 0$ и ${{d}_{k}} > 0$ от $m$ не зависят). Получаем теперь Значит и условие 2) тоже выполнено. Следствие установлено.

Замечание 1. Смысл свойства ненасыщаемости численного метода состоит как раз в следующем: если компакт ${{X}^{{{{k}_{2}}}}}$ устроен существенно проще, чем компакт ${{X}^{{{{k}_{1}}}}}$, то при прочих равных условиях вычислительный метод конструирует решение $\chi \in {{X}^{{{{k}_{2}}}}}$ уравнения (2.1) с точностью большей, чем в случае, когда $\chi \in {{X}^{{{{k}_{1}}}}}$.

Дискретизация (3.4) задачи (2.1) не создает коллизий точности вычислений с объемом перерабатываемой битовой информации, если решение $\chi $ достаточно гладкое, в частности, бесконечно дифференцируемое. При этом не возникает проблем, связанных с решением линейных алгебраических систем “большой” размерности. А утверждение теоремы 1 о том, что число обусловленности $\varkappa (I + A)$ не зависит от $m$, обеспечивает устойчивость итерационного процесса решения задачи (3.4). Иными словами, ошибки округления отыскиваемого решения $\bar {\chi }$ имеют примерно тот же порядок малости, что и ошибки, допущенные при вычислении матрицы $I + A$ и правой части $F$ (см. оценку (3.5')). Указанный факт является важным преимуществом итерационного метода решения линейной системы уравнений (3.4): хотя матрица $I + A$ не имеет специальной структуры и полностью заполнена, тем не менее существует эффективный итерационный метод ее численного решения.

Теорема 2. Последовательность функций

получаемая при решении системы (3.4) методом итераций, сходится к решению $\bar {\chi }$ задачи (3.4) так же, как геометрическая прогрессия со знаменателем, меньшим единицы.

Для окончательного заключения о близости приближения ${{\bar {\chi }}^{k}}$ к решению $\bar {\chi }$ задачи (3.4) требуется эффективный критерий прерывания итерационного процесса (3.6). Естественно принять за такой критерий количество верных значащих цифр в конструируемом ответе. Это означает, что, исходя из оценки ${{\left| {{{{\bar {\chi }}}^{k}} - \bar {\chi }} \right|}_{\infty }}{\text{/}}{{\left| {\bar {\chi }} \right|}_{\infty }} \leqslant \varepsilon $ при заданном $\varepsilon > 0$, проводить итерации следует до выполнения этого неравенства. Но априори само решение $\bar {\chi }$ нам неизвестно и указанный критерий рассматривать как эффективный нельзя. Укажем другой критерий, использующий последовательность значений ${{\bar {\chi }}^{k}}$ и число обусловленности $\varkappa (I + A)$, характеризующий трудность [5] решения задачи (3.4).

Теорема 3. Если

то

Доказательство этой теоремы простое и мы его приведем. Пусть $B = I + A$, тогда $B\bar {\psi } = F$ и, следовательно,

Учитывая, что получим Далее Теорема 3 доказана.

Замечание 2. Осуществить итерационный процесс (3.6) с любой заданной точностью $\varepsilon > 0$, используя только компьютерную конечно-разрядную арифметику, нельзя, если не предпринимать специальных мер [5]. Однако в тех случаях, когда число $\varkappa (I + A)$ невелико, получаются достоверные данные о величине погрешности.

Преимущества предлагаемого в теоремах 1–3 метода численного решения интегрального уравнения (2.1) сохраняются лишь при условии, что точность приближенной реализации интегрального оператора (2.2) имеет тот же порядок, что и величина ${{\left| \varrho \right|}_{\infty }} = {{\left| {JK[\chi - {{Q}_{m}}\chi ]} \right|}_{\infty }}$ при $m \geqslant {{m}_{0}}$ (см. (3.3)). И потому вопрос о ценности, или, что то же самое, о точности конструируемого компьютерного ответа, всецело находясь во власти способа численной реализации интегрального оператора $K$, напрямую зависит от качеств используемых квадратурных формул.

Однако сам вид выражения (2.2) показывает, что далеко не каждая квадратурная формула для вычисления элементов ${{a}_{{ik}}} = K[{{w}_{k}}]({{s}_{i}})$ представляет практический интерес. С одной стороны, ядро $k(s,\sigma )$ имеет на диагонали $\sigma = s$ “подвижную” логарифмическую особенность, обусловленную свойством модуля $q(\sigma ,s)$ эллиптических интегралов $E(q)$ и $D(q)$. С другой стороны, в точках $s$ вблизи оси симметрии $z$, являющейся зоной интенсивного роста функции $h_{ * }^{{ - 1}}(\sigma ,s)$, возникают явления пограничного слоя [13]. В полюсах же $\gamma (0)$ и $\gamma (\pi )$ выражение (2.2) вполне регулярно.

Требования к точности вычислений реализуются при использовании так называемых ненасыщаемых квадратурных формул. Такие формулы, учитывающие еще и специфику осесимметричной задачи (2.1), были построены в работе [14]:

Здесь ${{\tau }_{k}} = cos(\pi (2k - 1){\text{/}}2n)$, коэффициенты ${{c}_{k}}$ и ${{d}_{k}}$ указаны в [14], а функционалы погрешностей $\wp _{n}^{{C,D}}(f)$ квадратурных формул (3.7) оценены через характеристики ${{E}_{n}}(f)$ гладкости подынтегральной функции $f$ таким образом: $\left| {\wp _{n}^{C}(f)} \right|,\left| {\wp _{n}^{D}(f)} \right| \leqslant 4{{E}_{n}}(f)$; ${{E}_{n}}(f) = \mathop {inf}\limits_{{{P}_{n}} \in {{\mathcal{P}}^{n}}} \left\| {f - {{P}_{n}}} \right\|$ ($n > 0$ – целое) – наилучшее (чебышёвское) приближение непрерывной функции $f$ на отрезке $I \equiv [ - 1,1]$ многочленами из подпространства ${{\mathcal{P}}^{n}}$ алгебраических многочленов степени не выше $n - 1$.

4. ФЕНОМЕН НЕНАСЫЩАЕМОСТИ ВЫЧИСЛИТЕЛЬНЫХ ПРОЦЕДУР

Хотя математический смысл результатов, полученных в теоремах 1–3, прозрачен и прост, мотивация их практической применимости в цифровых вычислениях отнюдь не очевидна. Укажем доводы в пользу учета экстраординарных запасов гладкости решения $\chi (s)$ при приближенном решении уравнения (2.1). При этом оценим не только преимущества нового подхода, но и выявим аналитическую природу адаптивности процедур аппроксимации и квадратурных формул к запасам гладкости решения $\chi (s)$.

Пусть функция $\psi $ принадлежит $C_{ + }^{k}$, тогда ${{e}_{m}}(\psi ) \leqslant \tfrac{\pi }{2}\tfrac{{\left\| {{{\psi }^{{(k)}}}} \right\|}}{{{{m}^{k}}}}$ (теорема Джексона [2]). Пусть далее $\psi $ принадлежит гладкой шкале пространств $\bigcup\nolimits_{k \geqslant 0}^{} {C_{ + }^{k}} $ и ${\text{\{ }}G(k){\text{\} }}_{{k = 0}}^{\infty }$ – последовательность положительных чисел. С использованием условий

поставим в соответствие последовательности ${\text{\{ }}G(k){\text{\} }}_{{k = 0}}^{\infty }$ пару функций числового аргумента $x \geqslant 0$: Указанные классы ${{{\text{C}}}^{\infty }}$-гладких функций не пусты: им принадлежат, например, известные классы Жеврея, имеющие мажоранту $G(k) = {{A}^{k}}{{k}^{{\beta k}}}$ ($\beta \geqslant 1$, $A > 1$ – константы). С использованием этих определений теорема Джексона записывается в виде

Функция $\vartheta (m)$ в (4.1) является выражением нашего интуитивного представления о порядке убывания к нулю с ростом параметра $m$ аппроксимационных характеристик ${{e}_{m}}(\psi )$, а величина $\mu (m)$ при этом оказывается точностью приближения функции $\psi $.

Теорема 4 (см. [4]). При $x \geqslant 1$ функция $\vartheta (x)$ целочисленна, неотрицательна, не убывает, непрерывна справа и стремится к бесконечности при $x \to \infty $. Функция $\mu (x)$ строго монотонно убывает, всюду непрерывна и стремится к нулю при $x \to \infty $.

Следствие 2. При $p \geqslant 0$ верно следующее равенство: $\mathop {lim}\limits_{x \to \infty } {{x}^{p}}\mu (x) = 0$.

Из оценки (3.5), в силу неравенства (4.1) и теоремы 4, следует, что с ростом запаса гладкости решения $\chi $ задачи (2.1) скорость убывания к нулю погрешности построенного ненасыщаемого численного метода только возрастает. Причем с ростом $m$ метод самосовершенствуется и экспоненциальная его сходимость, в силу следствия 2, достигается на классе $C_{ + }^{\infty }$-гладких решений $\chi $. Так что, если $\chi $ из $C_{ + }^{\infty }$ и $G(k) = {{A}^{k}}{{k}^{{\beta k}}}$ при $A > 1$ и $\beta \geqslant 1$, то ${{e}_{m}}(\chi ) \leqslant c{{e}^{{ - r\sqrt[\beta ]{m}}}}$, где $c$, $r$ – положительные константы.

Попытка понять природу приспособляемости квадратурных формул (3.7) к запасам гладкости подынтегральных функций $f$ привела, в свою очередь, к некой новой формализации оценок их функционалов погрешностей с помощью специально сконструированных для этого монотонных функций $\lambda (x)$ и $\theta (x)$ числового аргумента $x$.

Итак, пусть $C[I]$ – пространство непрерывных на отрезке $I \equiv [ - 1,1]$ функций с чебышёвской нормой $\left\| {\, \cdot \,} \right\|$, а ${{C}^{k}}[I]$ ($k \geqslant 0$ – целое число) – пространство $k$ раз непрерывно дифференцируемых на $I$ функций. При этом для функции $f$ из ${{C}^{k}}[I]$ справедлива теорема Джексона–Зинвела [15]:

Пусть

Определим следующие функции числового аргумента $x \in [0,\infty )$: В этих обозначениях неравенство Джексона–Зинвела примет следующий вид:

Теорема 5 (см. [14]). При $x \geqslant 1$ функция $\theta (x)$ целочисленна и неотрицательна, не убывает, непрерывна справа и стремится к бесконечности вместе с $x$. Функция $\lambda (x)$ строго монотонно убывает, непрерывна справа и стремится к нулю при $x \to \infty $. Функция $\lambda (x)$ имеет разрывы слева лишь в точках разрыва функции $\theta (x)$.

Следствие 3. При $p \geqslant 0$ верно предельное равенство $\mathop {lim}\limits_{x \to \infty } {{x}^{p}}\lambda (x) = 0$.

Теорема 5 расширяет сферу практической применимости полиномиальной аппроксимации ${{C}^{\infty }}$-гладких функций. Действительно, рассматривая ${{C}^{\infty }}$-гладкие на отрезке $I \equiv [ - 1,1]$ функции, мы считали их бесконечно дифференцируемыми всюду вплоть до границ $I$. Представляет интерес и тот случай, когда производные функции $f$, будучи ограниченными в замкнутом промежутке $I$, допускают вблизи его концов рост, характер которого фиксирован с помощью некоторого числового параметра. Подобные ситуации типичны для вычислительной практики. Они возникают, например, при численном интегрировании функций, имеющих в режимах своего поведения резкие переходные зоны – “пограничные слои”. Размеры этих переходных зон обычно характеризуются величиной числового параметра ${{\varepsilon }_{0}}$, $0 < {{\varepsilon }_{0}} < 1$, называемого толщиной пограничного слоя. Может ли это (если может, то как) обязательно отразиться на сходимости к нулю с ростом числа узлов $n$ функционалов $\wp _{n}^{{C,D}}(f)$ погрешностей квадратурных формул (3.7)?

Проведенные К.И. Бабенко изыскания [16] сформировали представление о том, в каких именно конструктивных терминах удобно описывать указанную специфику задач и за счет каких ресурсов возможна ее компьютерная численная нейтрализация. Дадим более четкую формулировку этой проблемы, введя специальный термин, характеризующий рост градиентов функций вблизи концов отрезка $I$.

Определение (см. [14]). Функция $f(\tau )$, принадлежащая ${{C}^{\infty }}[I]$, имеет на $I \equiv [ - 1,1]$ пограничный слой толщиной $0 < {{\varepsilon }_{0}} < 1$, если существуют малое положительное число $\eta = \eta ({{\varepsilon }_{0}})$ и положительная функция $F(k)$, не зависящая от ${{\varepsilon }_{0}}$, такие что при любом целом $k \geqslant 0$ справедливо неравенство

Наводящим соображением к выяснению характера влияния пограничного слоя на величины погрешностей $\left| {\wp _{n}^{{C,D}}(f)} \right|$ квадратурных формул (3.7) будет служить следующий [17], восходящий к работе С.М. Никольского [18], результат.

Теорема (В.К. Дзядык). Для того, чтобы $k$-я производная функции $f$ удовлетворяла на отрезке $I \equiv [ - 1,1]$ условию Гëльдера с показателем $0 < \gamma < 1$, необходимо и достаточно, чтобы при любом целом $n \geqslant k$ существовал такой алгебраический многочлен ${{P}_{n}} \in {{\mathcal{P}}^{n}}$ степени $n - 1$, что для всех $t \in I$ справедливы оценки

где ${{A}_{k}}$ – постоянная, не зависящая от $t$ и $n$.

Теорема Дзядыка усиливает неравенство Джексона–Зинвела (4.2), ибо при равномерной оценке ${{E}_{n}}(f) \leqslant M{{n}^{{ - (k + \gamma )}}}$, устанавливает возможность приближения функции $f \in {{C}^{{k + \gamma }}}[I]$ вблизи концов отрезка $I$ c погрешностью $O({{n}^{{ - 2(k + \gamma )}}})$.

Следствием указанной фундаментальности, существующей в природе полиномиальной аппроксимации гладких функций на отрезке $I$, служит следующая

Теорема 6 (см. [14]). Если $f$ принадлежит ${{C}^{\infty }}[I]$ и выполнено (4.3) и $\eta = {{\varepsilon }_{0}}{\text{/}}2$, то

Коэффициенты ${{F}_{k}}$ вычисляются по заданным значениям $F(k)$ из (4.3).

Прикладное значение оценки (4.4) для ненасыщаемых квадратурных формул (3.7) заключается в следующем: за счет перераспределения пограничного слоя по всему отрезку $I$ его толщину удается увеличить до значения $\sqrt {{{\varepsilon }_{0}}} $. Этот факт и лежит в основе способа численной нейтрализации пограничного слоя толщиной ${{\varepsilon }_{0}} > 0$.

В самом деле, согласно теореме 5, выполнимость неравенств

представляющих условия нейтрализации пограничного слоя толщиной ${{\varepsilon }_{0}}$, обеспечивается выбором в формулах (3.7) числа узлов $n$, превышающих некоторое пороговое значение $n > {{n}_{{min}}}({{\varepsilon }_{0}})$. Среди оценок (4.5), отвечающих различным $0 \leqslant k \leqslant n$, имеется наилучшая, номер которой ${{k}_{0}} = \theta (n)$ – порядок максимальной производной, содержащейся в этой оценке; производные порядков $k > {{k}_{0}}$ могут влиять на величину самой оценки лишь в случае $n > {{n}_{{min}}}$, поэтому нейтрализация пограничного слоя квадратурными формулами (3.7) осуществляется за счет выбора числа узлов $n > {{n}_{{min}}}$.

Следовательно, наличие у подынтегральной функции больших запасов гладкости создает благоприятные предпосылки для эффективного построения числового ответа, причем интерес к ${{C}^{\infty }}$-гладким функциям вовсе не кажется неестественным [6].

Квадратурные формулы с главным членом погрешности (т.е. насыщаемые – формулы прямоугольников, трапеций, Симпсона и ряд других: см. [2]) способностью к нейтрализации пограничного слоя не обладают; для них пограничный слой является “камнем преткновения” при их компьютерной реализации.

Аппроксимация интегрального оператора $K$ сводится к применению формул (3.7). Покажем, как это сделать. Представление (2.2) преобразуется к виду (см. [13])

в котором логарифмическая особенность выделяется явно, а $A(\sigma ,s)$, $B(\sigma ,s)$ – достаточно гладкие и равномерно непрерывные в области $[0,\pi ] \times [0,\pi ]$ функции.

Зафиксировав в (4.6) параметр $s \in (0,\pi )$ и заменив неявно переменную $\sigma $ новой переменной $\tau \equiv \tau (\sigma ,s) = {{sin\tfrac{{(\sigma - s)}}{2}} \mathord{\left/ {\vphantom {{sin\tfrac{{(\sigma - s)}}{2}} {sin\tfrac{{(\sigma + s)}}{2}}}} \right. \kern-0em} {sin\tfrac{{(\sigma + s)}}{2}}}$, получим

Здесь подвижная логарифмическая особенность в (4.6) перешла в неподвижную – середину отрезка $I$. А для функций $A(\sigma ,s)$ и $B(\sigma ,s)$, равномерно непрерывных на квадрате $[0,\pi ] \times [0,\pi ]$, выполняются при $k \geqslant 0$ соотношения: где $\tilde {g} \equiv \tilde {g}(\tau ) = g(\sigma (\tau ,s),s)$, а в качестве $g(\sigma ,s)$ используются функции $A(\sigma ,s)$ или $B(\sigma ,s)$; функция $\sigma (\tau ,s)$ является обратной к функции $\tau (\sigma ,s)$ и $\Psi (\tau ) = \chi (\sigma (\tau ,s))$.

Пограничный слой в (4.7) толщиной ${{\varepsilon }_{0}} = 0.5sins$ выделен, в силу (4.8), явно. И потому вычисление элементов ${{a}_{{ik}}} = K[{{w}_{k}}]({{s}_{i}})$ матрицы $A$ квадратурными формулами (3.7) с числом узлов $n > {{n}_{{min}}}({{\varepsilon }_{0}}) > m \geqslant {{m}_{0}}$ позволяет его, согласно неравенству (4.5), эффективно нейтрализовать. При этом матрица $A = ({{a}_{{ij}}})$ вычисляется с любой заданной точностью $(1 + 2\left\| {{{Q}_{m}}} \right\|){{e}_{m}}(\chi )$, поскольку из теоремы 4 (в силу оценки (4.1)) следует, что с ростом запаса гладкости решения $\chi $ задачи (2.1) скорость убывания к нулю характеристик ${{e}_{m}}(\chi )$ с ростом $m$ только возрастает. Поэтому в отличие от методов, имеющих главный член погрешности, построенный метод с ростом параметра $m$, согласно теореме 4, самосовершенствуется и, преодолев барьер степенной сходимости, максимума своей эффективности – экспоненциальной сходимости (или сверхсходимости) – достигает, согласно следствию 2, на классе $C_{ + }^{\infty }$-гладких решений $\chi $ задачи (2.1). Так, если $\chi \in C_{ + }^{\infty }$ и $G(k) = {{A}^{k}}{{k}^{{\beta k}}}$ при $A > 1$ и $\beta \geqslant 1$, то ${{e}_{m}}(\chi ) \leqslant c{{e}^{{ - r\sqrt[\beta ]{m}}}}$, где $c$, $r$ – положительные константы.

Таким образом, информация о бесконечной гладкости решения $\chi $ задачи (2.1) обретает на практике особую важность, что принципиально отличает ненасыщаемый численный метод от метода с главным членом погрешности, т.е. насыщаемого. Сверхсходимость ненасыщаемого метода решения уравнения (2.1) обеспечивает экономность конструирования компьютерного числового ответа с точностью, определяемой исключительно аппроксимативными возможностями компакта ${{C}^{\infty }}$-гладких функций, содержащего точное решение $\chi $ уравнения (2.1).

5. ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ СПЕКТРАЛЬНОЙ ЗАДАЧИ ДЛЯ УРАВНЕНИЯ (2.1)

В отборе итерационного метода решения системы уравнений (3.4) следует быть весьма осмотрительным, с тем чтобы не лишить метод практической целесообразности – допустимого по точности ответа. Одним из факторов, обеспечивающих эффективное его функционирование, является удачный выбор итерационных параметров. Причем сам выбор не является волюнтаристским: содержательность его всякий раз должна подкрепляться информацией, “скрывающейся” в спектральном портрете матрицы $A$ и всецело определяемый способом дискретизации задачи (2.1): наличие у $A$ кратных или достаточно близких собственных значений может серьезно сказаться на скорости сходимости итерационного процесса, ибо поведение степеней ${{A}^{k}}$ определяется структурой спектра матрицы $A$. И худший случай, который может здесь представиться – это наличие у матрицы $A$ жордановых клеток (отсутствие у $A$ этих патологий отнюдь не очевидно, ввиду несамосопряженности задачи (2.1)).

Эту трудность удается преодолеть, ориентируясь на ключевые свойства оператора $K$ задачи (2.1). Оператор $K$ – компактен, ненулевые точки его спектра, т.е. вещественные простые полюсы резольвенты $R(\zeta ,K) \equiv {{(K - \zeta I)}^{{ - 1}}}$, содержатся в промежутке $\left( { - 1,1} \right]$, а его собственные функции $\psi $ непрерывны [8].

Исходя из этого, выявим субстрат тех общих представлений, который позволяет автоматически, исходя из конкретно складывающейся ситуации, отбирать параметры итерационного процесса (3.6). Иными словами, покажем, что ненасыщаемая дискретизация задачи (2.1) приводит к алгебраической задаче (3.4) с “хорошей” матрицей $A$, наследующей спектральные свойства оператора $K$. И потому в ее спектральном портрете указанные патологии невозможны, если они отсутствуют у оператора $K$.

Действительно, рассмотрим спектральную задачу для оператора $K:{{C}_{ + }} \to {{C}_{ + }}$:

Здесь число $\lambda $ – это в точности искомое собственное значение задачи (5.1), а $\psi $ – соответствующая числу $\lambda $ собственная функция $\psi (s)$.

Далее, осуществляя по известной схеме дискретизацию задачи (5.1), получаем

Здесь $\nu $ и $\bar {\psi } \in {{\mathbb{R}}^{{m + 1}}}$ – соответственно собственное значение и собственный вектор матрицы $A:{{\mathbb{R}}^{{m + 1}}} \to {{\mathbb{R}}^{{m + 1}}}$, а компоненты $\bar {\psi }$ – приближенные значения в узлах ${{s}_{j}}$, $0 \leqslant j \leqslant m$, собственной функции $\psi $ задачи (5.1) (матрица – $A$ та же, что и в (3.4)).

Оценим возмущение, вносимое в собственное число $\lambda $ спектральной задачи (5.1) отбрасываемой погрешностью – вектором $\delta = - JK[\psi - {{Q}_{m}}\psi ]$, $\delta \in {{\mathbb{R}}^{{m + 1}}}$.

Идея решения этого вопроса заимствована нами из работы [19]. А реализовать ее удалось, используя теорию регулярных возмущений линейных ограниченных операторов в банаховом пространстве [20]. Согласно этой теории, возмущение спектра оператора $K$ связывается с равномерной (по норме) сходимостью последовательности приближающих $K$ операторов $K{{Q}_{m}}$. Выбор интерполяционного проектора ${{Q}_{m}}$ в форме (3.2) обеспечивает выполнение нужного свойства. Справедлива

Теорема 7 (см. [7]). Последовательность операторов ${\text{\{ }}K{{Q}_{m}}{\text{\} }}$ сходится при $m \to \infty $ равномерно к компактному оператору $K$ задачи (5.1):

Указанные спектральные свойства оператора $K$ задачи (5.1) в сочетании с введенным С.Л. Соболевым [21], [22] понятием близости операторов, отличным от близости по норме, дают возможность оценить, насколько близки спектры $\sigma (K)$ и $\sigma (K{{Q}_{m}})$ двух несамосопряженных операторов, если близки сами операторы $K$ и $K{{Q}_{m}}$.

Именно, выясним, как изменяются собственное значение $\lambda $ и собственная функция $\psi $ оператора $K$ задачи (5.1) при малых (в равномерной норме) его возмущениях.

Теорема 8 (см. [7]). Пусть внутри гладкого замкнутого контура ${{\Gamma }_{\lambda }}$ на комплексной плоскости находится ровно одно собственное значение $\lambda $ задачи (5.1). Тогда внутри ${{\Gamma }_{\lambda }}$ имеется ровно одно собственное число $\nu $ матрицы $A$ задачи (5.2). Справедлива также оценка погрешности

и положительная постоянная ${{c}_{4}}$ вычисляется эффективно.

Следствие 4. Матрица $A$ имеет собственный вектор $\bar {\psi }$ такой, что

и положительная постоянная ${{c}_{5}}$ эффективно вычисляется.

Таким образом, число $\nu $ и полином ${{Q}_{m}}(s;\bar {\psi })$ – это искомые приближения к первому собственному числу $\lambda $ и собственной функции $\psi $ спектральной задачи (5.1).

Из теоремы 8 следует, что скорость сходимости последовательности приближенных решений $\nu \equiv {{\nu }_{m}}$ и $\bar {\psi } \equiv {{\bar {\psi }}_{m}}$ спектральной задачи (5.2) с ростом параметра $m \geqslant {{m}_{0}}$ определяется исключительно гладкостью собственной функции $\psi $ задачи (5.1). Причем возмущение, вносимое построенной ненасыщаемой дискретизацией в спектральную проблему (5.1), зависит от близости резольвент операторов $K$ и $K{{Q}_{m}}$ вблизи первого собственного числа оператора $K$. При этом, если простое собственное значение $\lambda $ оператора $K$ хорошо отделено от остальных его собственных значений, то построенный нами итерационный процесс (3.6) численно устойчив.

Замечание 3. Элементы матрицы $A = \left( {{{a}_{{ik}}}} \right)$, ${{a}_{{ik}}} = K[{{w}_{k}}]({{s}_{i}})$, $0 \leqslant k \leqslant m$, $0 \leqslant i \leqslant m$, спектральной задачи (5.2) всегда определены с точностью, которая не может быть меньше ${{\left| \delta \right|}_{\infty }}$. Поэтому мы можем изменять элементы ${{a}_{{ik}}}$ матрицы $A$ на величину порядка ${{\left| \delta \right|}_{\infty }} = {{\left| {JK[\psi - {{Q}_{m}}\psi ]} \right|}_{\infty }}$. Иначе говоря, задача абсолютно точного определения собственных чисел матрицы $A$ лишена всякого смысла [23]. В общей ситуации имеет смысл лишь задача о “почти собственных значениях” матрицы $A$ при величине невязки порядка ${{\left| \delta \right|}_{\infty }}$. Любой метод численного решения алгебраической задачи на собственные значения есть метод вычисления “почти собственных значений” и потому важно уметь адаптировать его к классу корректности рассматриваемой задачи. Из сказанного явствует, насколько важно осуществлять ненасыщаемые дискретизации спектральных задач.

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

Чтобы оценить максимум модуля тригонометрического многочлена $T(s)$ степени $m$, вычислим сначала его значения в равноотстоящих узлах ${{s}_{j}} = \tfrac{{2\pi j}}{{2m + 1}}$, $j = 0,\; \ldots ,\;2m$. Затем, воспользовавшись неравенством Лебега [2], запишем оценку

Далее, по теореме С.Н. Бернштейна [2], имеем $\left| {{{T}^{{(k)}}}} \right| \leqslant {{m}^{k}}H$.

Пусть $M \gg m$. Вычислим значение $T(s)$ по бóльшему числу узлов ${{\sigma }_{j}} = \tfrac{{2\pi j}}{{2M + 1}}$, $j = 0,\; \ldots ,\;2M$. Если ${{\sigma }_{ * }}$ – точка максимума многочлена $\left| {T(s)} \right|$, ${{\sigma }_{l}} \leqslant {{\sigma }_{ * }} \leqslant {{\sigma }_{{l + 1}}}$, то по формуле Тейлора имеем в некоторой точке $\sigma $, ${{\sigma }_{l}} \leqslant \sigma \leqslant {{\sigma }_{ * }}$:

Выбирая число $M$ подходящим образом, можем оценить $\mathop {max}\limits_{0 \leqslant s < 2\pi } \left| {T(s)} \right|$ по узловым значениям $\left| {T({{s}_{j}})} \right|$, $j = 0,\; \ldots ,\;2M$. Учитывая, что в величине $\mathop {max}\limits_{0 \leqslant s < 2\pi } \left| {T(s)} \right|$ нужно лишь небольшое число значащих цифр, вычисления такого рода легко осуществить для полиномов $T(s)$ не очень больших степеней.

Методика численного решения задачи (2.1) протестирована нами на эллипсоиде вращения с удлинением, равным $1000$ (“игла”), как наиболее “трудном” для компьютерных расчетов примером. Задача ставилась таким образом, чтобы было получено большое количество верных десятичных разрядов в отыскиваемом решении задачи (2.1) при сравнительно небольших размерах матриц линейных алгебраических систем, к решению которых сводится задача. Численное решение задачи с 8–10 верными десятичными разрядами было получено при следующих, принятых в алгоритме, значениях числовых параметров: $m = 20$ (теорема 1), $n = 501$ (теорема 6) (см. [13], [24]).

Проведенные компьютерные расчеты свидетельствуют о преимуществе ненасыщаемых численных методов при решении эллиптических краевых задач, в то время как удлинение эллипсоида, равное 25, оказывается непреодолимым препятствием для любых насыщаемых (т.е. с главным членом погрешности) численных методов [2].

Таким образом, построенная в работе ненасыщаемая методика численного решения уравнения (2.1) с компактным оператором $K$ снабжает нас элегантным вычислительным средством, способным в дискретизованной форме наследовать как дифференциальные, так и спектральные характеристики оператора $K$. Это служит основанием для построения компьютерного числового ответа гарантированной точности [5], если решение $\chi $ уравнения (2.1) достаточно гладкое, например, ${{C}^{\infty }}$-гладкое.

Не стоит, однако, заблуждаться, полагая, что сколько-нибудь существенный прогресс для гладких трехмерных областей возможен простым переносом его с осесимметричного случая. Общий трехмерный случай все еще далек от теоретического изыска осесимметричного и требует привлечения глубоких фактов многомерного анализа: проблема ненасыщаемой аппроксимации функций на гладких многообразиях, гомеоморфных двумерной сфере, до сих пор эффективно не разрешена [25], [26]. И ее разрешение, как и предсказывает теория [1], [2], относится к сфере в уже гораздо бóльшей степени интеллектуальной, и на первый план выдвигается проблема конструирования ненасыщаемых кубатурных формул [27].