К проблеме численной реализации интегральных операторов осесимметричных краевых задач (алгоритмы без насыщения) Текст научной статьи по специальности «Математика»

ISSN 2074-1863 Уфимский математический журнал. Том 4. № 4 (2012). С. 22-37.

УДК 519.642+532.582.33

К ПРОБЛЕМЕ ЧИСЛЕННОЙ РЕАЛИЗАЦИИ ИНТЕГРАЛЬНЫХ ОПЕРАТОРОВ ОСЕСИММЕТРИЧНЫХ КРАЕВЫХ ЗАДАЧ (АЛГОРИТМЫ БЕЗ НАСЫЩЕНИЯ)

В.Н. БЕЛЫХ

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

Метод оттестирован на задаче прецизионного вычисления интеграла Гаусса теории гармонического потенциала в эллипсоидах вращения большого удлинения.

Ключевые слова: ненасыщаемый численный метод (unsaturated numerical method), интеграл Гаусса (Gauss integral), осесимметричная область (axisymmetric region), нена-сыщаемая квадратурная формула (quadrature formula without saturation).

1. Преамбула.

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

С расширением вычислительной практики и появлением в научном обороте новых классов прикладных задач [1] всё активнее начинает переосмысливаться роль запаса гладкости их решений, понимаемого как наличие у решений набора непрерывных производных, превышающего необходимый и вполне определенный порог. Но такие экстраординарные запасы гладкости обычно потенциальны и не реализуются, если у погрешности метода имеется главный член, так как с расширением аппроксимационных возможностей компакта решений, определяемых гладкостью его элементов, погрешность метода не уменьшается [2]. Указанный дефект сильно снижает практическую ценность численных методов, имеющих главный член погрешности. Он же и объясняет тенденцию, почему вычислители обычно стремятся использовать методы высоких порядков, в которых отмеченный дефект, хотя и не исчезает полностью, но уменьшается.

Современная вычислительная наука (как фундаментальная, так и прикладная) испытывает насущную потребность в появлении численных методов, лишённых указанного

V.N. Belykh, The problem of numerical realization of integral operators of axisymmetric boundary value problems (algorithms without saturation). © Белых В.Н. 2012.

Работа выполнена при финансовой поддержке РФФИ грант 11-01-00147-a и Междисциплинарного интеграционного проекта № 40 Президиума СО РАН. Поступила 22 сентября 2011 г.

дефекта - наличия главного члена погрешности, - и названных К.И. Бабенко ненасыща-емыми [2]. Генетически эти методы примыкают к методам переменных порядков, отличаясь от них, лишь повышенной согласованностью с классом корректности задачи: с ростом запаса гладкости решений ненасыщаемые численные методы самосовершенствуются, т.е. автоматически черпают приращение своей практической эффективности непосредственно в дифференциальной природе отыскиваемых решений задач (феномен ненасыщаемости [3]). При этом пик эффективности - экспоненциальная сходимость - достигается на классах С -гладких решений.

Существуют классы задач, в которых информацией об экстраординарной гладкости решений можно распорядиться с пользой для вычислений. Так, в эллиптических задачах качественные дискретизации отличает их способность к наследованию характеристического свойства - шаудеровских оценок гладкости [4]. Ранее эта способность едва ли была востребована практикой, но актуальной её делает известный контекст - стремление к числовому ответу гарантированного качества. Это в совокупности с хорошей обусловленностью дискретизации способно приводить к истинной terra incognita компьютерных вычислений -их доказательности [2,5]. Хорошая обусловленность дискретизации при этом совершенно необходима, поскольку математические сущности, пребывающие на этапе создания численного метода, очень часто становятся идеальным полем для всевозможных спекуляций, не принимающих во внимание реальности компьютерных вычислений - ошибок округления.

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

Между тем использование стандартных, т.е. с главным членом погрешности, кубатур-ных формул сталкивается лишь с формальной констатацией существующих трудностей ("проклятий") финитной математики (размерности пространства независимых переменных [2,6] и кратности интеграла [7]) без указания возможных путей их преодоления. В результате приходится надеяться только на удачу, а не на гарантированный успех. Именно по этой причине проблема отыскания гарантированного качества численного решения уравнения Лапласа в гладких трёхмерных областях достаточно произвольной формы до сих пор представляется очень трудной вычислительной проблемой. Разрешение её, как и предсказывает теория [2], относится к сфере в уже гораздо большей степени интеллектуальной: на первый план выдвигается проблема построения ненасыщаемых кубатурных формул. Существующие на этом пути трудности усугубляются и отсутствием тех прозрений здравого смысла, которые бы могли удовлетворить реальным запросам практики. Впрочем, в недавней работе [8] была предпринята интересная попытка разобраться в этом.

Высказанная точка зрения на проблему качества компьютерных вычислений оказала серьёзное влияние и на выбор материала для статьи. В ней мы ограничимся рассмотрением осесимметричных краевых задач, предъявив новую - ненасыщаемую - методику дискретизации интегральных операторов этих задач. В общих чертах она описана в [5]. Здесь же на примере прецизионного вычисления интеграла Гаусса обрисован общий источник вычислительных трудностей этих задач: сильный рост интегрантов вблизи оси

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

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

2. Интеграл Гаусса на осесимметричной поверхности

Новые численные методы создаются под влиянием или с целью разрешения какой-то вполне конкретной прикладной задачи. В исследуемой нами проблематике эту роль полноценно исполняет интеграл Гаусса теории гармонического потенциала [9]. Значение этого интеграла не зависит от поверхности, по которой он вычисляется. В то же время его аналитическая природа достаточно сложна, чтобы служить моделью для общей задачи в гладкой осесимметричной области. И в этом - его уникальность, как тестового примера.

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

Пусть х € R3, х = (х,у, z); ш С R3 - область с осью симметрии z, ограниченная С ^-гладкой замкнутой поверхностью вращения дш; величины г = л/х2 + у2 их- инварианты группы вращения области ш относительно оси z. Меридиональное сечение дш -это параметризованная кривая 7 : [0,1] ^ [г( s), z(s)}, г > 0, dz/ds > 0, 7(5) € Сте[0,1]. Точки 7(0) и 7(1) - это полюсы дш, а r(s) и z(s) имеют С ^-гладкие 2-периодические продолжения с [0,1] на [0, 2] (нечётное, чётное соответственно).

Кривую 7(s) будем считать выпуклой, так что точки, близкие к оси z, всегда сосредоточены вблизи полюсов 7(0) и 7(1) поверхности дш. Общий (невыпуклый) случай исследуется соответствующим разбиением ( ) на части.

Положение точек х и ^ на дш определим цилиндрическими координатами (r,z, v) и (р,(,р) соответственно: х = x(s, v) = (reos v, rsin v, z), £ = x(a,p) = (pcosp, psinp, (). Здесь обозначено: г = r(s), z = z(s), p = r(a), ( = z(a) и 0 < s,a < 1, 0 < v,p < 2ж.

Условимся считать точку х € дш фиксированной, а точку ^ € дш - переменной.

Нормаль п и элемент площади dш^ в точке ^ € дш обозначим соответственно:

п = ^ — (' cos р, — í-1 (' sin р, í-1 p'^j , dш% = pö da dtp,

где (' = d(/da, p' = dp/da, 5 = 5(a) = \Jp'2 + ('2.

Прямое значение нормальной производной ~Wx¡ < п > функции /(£, х) на поверхности вращения дш определим равенством:

1 Í ,,Bf ,df\ lC Тр

Vxi < П >= —8 1 (С' .

х е дш

Рассмотрим интеграл Гаусса по замкнутой поверхности вращения дш:

1 Г 1

ж = Г(ж) = - V< п > dш£, х е дш, £ е дш. (1)

дш

Здесь Р = | £ — х | - расстояние между точками х е дш, причём Р2 = г2 + р2 + (( — z)2 — 2 pr cos (ip — v), а интеграл понимается в смысле главного значения.

Введём обозначения:

h2 = h2(a, а) = (p - г)2 + (С - z)2, К = hl(a, s) = (p + r)2 + (( - z)2, q = q(a, s) = 4prh-2, u ^ 1 -q = h2 h-2,

7T/2 7T/2

X(/3) = J (1 -P sin2 0)"1/2 d0, E(P) = j (1 - P sin2 e)1/2 <№, D(P) = K(P) - E(P)■ 0 0 Здесь 0 < P < 1 - модуль полных эллиптических интегралов K(P), E(P), D(P).

Осевая симметрия поверхности дш позволяет один раз в (1) проинтегрировать, сведя вычисление в меридиональную плоскость v = const = 0 на основе следующих фактов. ЛЕММА 1. На гладкой замкнутой поверхности вращения дш верно представление:

„ 1 s h2 С' 1 - cos ю W— = Q(a, s)--г----,

V X р а) р з ' ^ р з )

где обозначено:

(a - s)-1 п • r[a, s] / | r[a, s]| 2, a = s, 0.5 (z V - r'z'') / d3, a = s,

t

Q = Q(a, s)- r „,. .,,w,3

r[a, s]

p - p - -cos <p, - sin <p.

d = d( s ) = V r'2 + z'2.

a — s a — s a — s

Здесь r', z' и r", z" - первые и вторые производные функций г = r(s), z = z(s) по локальной координате s соответственно.

ЛЕММА 2. На гладкой замкнутой поверхности вращения дш справедливо:

2ж 2ж

j Р-3 d^p = 4 h2 h-lE(q), j (1 — cos p)P-3 dip = 2 r-1 p-l h-1 D(q). 0 0

С учётом лемм 1,2 интеграл Гаусса (1) преобразуется к виду:

1

„ = Ц.)^^, s) E (q) — cm) КЧ,, [0, я

0

Наличие здесь полных эллиптических интегралов требует специальных методов аппроксимации E(q) и D(q) с учётом того, что модуль q = q(a, s) является функцией двух точек: переменной a и фиксированной s. При этом всегда q(s, s) = 1. В [10] получены представления

E(cd = Ep(Q) — е*Р ((1) lnu, d(q) = D^*(Q) — d*P (0) lnu.

Здесь параметр р > 0 - целое число, алгоритм выбора которого будет определён позже; и = 1 — q, а функции E*(q), e*(q), D*(q), d*(q) задаются формулами:

E *()=,0.5n(1 — qE о <q>) + e*p (q)ln(1 — q), если 0 <q< 0.5, v ' —Dp <и> up+l In и + 1 + uG < и >, если 0.5 < q< 1,

I

I

где

D *() = j'KqD 0 <q> +d * (q) ln(1 - q), если 0 <q< 0.5,

p(q) = ^0.5Ep <u>up+1 In и + \n(4/e) -uW<u>, если 0.5 < q< 1.

ж n—n— 1 P n ^

TP ^ ^ X^ ^nXП P 1 \ - n InXП ^

EP <X>= x^ 2n- 1 , 6p (x) = ¿^ 2n- 1 , e*(x) = 0, G<x>=^Gn,

n=p+1 n= 1 n= 1

те п—р—1 Р п

бр <х>= £ ' &» = -0-5£ ' &0*(х) = 1,

п=р+1 п=1 п=1

п 1пхп 1 / 1 \ 1 пхп 1

л п + л п-1) 1 , ^ + л п) п

2п - 1 ' ' \2п - 1 ) 2п - 1 '

7 о = 1, 7п =((1/2) п/п^ , Л о = 1п 4, Лга = ф(п + 1) - ф(п + 1/2), п> 0,

(а)п = Г(а + п)/Г(а), а> 0, п> 0, ф(х) = &Г(г)/&г (г > 0).

Здесь Г(г) и ф(г) соответственно Г- и ^-функции Эйлера (не путать с интегралом Г(з)). Алгоритмы вычисления конечных сумм е* (д) и & * (д) и быстросходящихся степенных рядов Е0 < д >, О0 < д >, Ер < и >, Ор < и >, Ш < и >, С < и > указаны в [10].

Таким образом интеграл Гаусса (1) приводится к стандартному для осесимметричных областей виду:

1

я- = Г(*) = У (^2р6 П(а, з) Е*р(д) - (' О *р(д)^ к-1 &а-

1а. (2)

У ^2 р5П(а, з ) е*р (д) - (' &*р (д)^ к-1 1п(1 - д)&а. о

Обратим внимание на то, что в полюсах 7(0) и 7(1) гладкой поверхности вращения дш представление (2) вполне регулярно, т.е. не содержит логарифмической особенности.

3. Специфика вычислительных трудностей в осесимметричных задачах

Далеко не каждая квадратурная формула для приближённой реализации (2) может представить практический интерес: структура сложности интегранта в (2) определяется как свойствами модуля д(а, 8), так и поведением функции к-1 (а, 8) в точках, близких к оси симметрии , т.е. вблизи полюсов (0) и (1). Действительно, на диагонали а = инте-грант в (2) имеет "подвижную" логарифмическую особенность, а в точках з, близких к оси симметрии г, обнаруживает ярко выраженный рост - пограничный слой [3,11]. Математически это означает наличие в (2) параметра при стремлении к нулю которого, выявляется структура роста интегранта и его производных в точках, близких к оси симметрии г. В (2) этим параметром служит весовая функция к-1(а, 8). Тем самым, явившись своеобразной платой за цилиндрическую симметрию задачи, пограничный слой в (2) оказывается присущ всем задачам с осевой симметрией. Понятно, что это не может не сказаться на реализации (2) стандартными квадратурными формулами.

Прецизионные вычислительные эксперименты показывают, что для эффективной компьютерной реализации представления (2) необходимо:

1) распорядиться свойствами модуля д(а, 5 ) для равномерного по (а, 5 ) выделения логарифмической особенности, и общепринятое её выделение по правилу

1п(1 - д) = 21п | а -в| + А(а, в), (а, в) е [0,1] х (0,1)

непригодно вблизи оси г, так как функция А(а, в) не является равномерно непрерывной;

2) обеспечить численную нейтрализацию пограничного слоя.

Математические изыскания [3] сформировали представление о том, в каких именно конструктивных терминах удобно описывать указанные трудности и за счёт каких ресурсов возможна их численная нейтрализация.

1

Преобразуем (2) к виду, удобному для дальнейшего анализа. Введём функции, равномерно непрерывные на квадрате К х К (здесь К = [0,1]):

R V, S) =( ■ Р—Is)) +( Л(as)) , R l(a, S) =( .P::+S)) +( Л(als)) ,

(r/ sin жs)(p/ sin na) R2{a, s)

Q(a, ¿0 = 4-^-г-, В(a, s) = —-г. (3)

R*2( Т, ) R*2( Т, )

Функции в (3) равномерно непрерывны и принадлежат С™[К х К] при 'y(s) Е С ™[К]. Интеграл (2), в котором переменная a неявно заменяется новой переменной

/ \ / \ • ж(а — s) , . ж (a + s) t = ts(a) = t(a, s) = sin У — J / sin У ^ }, (4)

где s E (0,1) - фиксировано, а функция a = as(t) = a(t, s) является обратной к t = ts(a) = t(a, s), имеет следующий вид

i i

eT(s) = J ГС (t)d t — J TL(t)ln\t\d t, e = e(s) = 0.5n sinns, (5)

-i -i

где подынтегральные функции Г с(t) и ^b(t) задаются формулами (см. формулу (2))

Гс (t) = (2р Ü (Е * — 1) Ц\ ) — ¿-1( '(D * — 1) d* )) Ъ, TL(t) = 2¡2p П Ц — Г1 С ' d* )) а,

a(t) = Ш-1 sin Ж (Т+ ^ , 6 = lnB. (6)

Для равномерно непрерывной на К х К функции g(a, s) по определению

д = g(t) = g(as(t), s), tEl = [ —1,1], sE (0,1) — фиксировано. (7)

Если g(a, s) E СГХ[К х К], то 7j(t) Е Си справедливо равенство

4~}g(t) = z-k{ sin2 (Т + S g(a, s), £ = £(s) = 0.5tt sinn s, Vk> 0. (8) dt J \ 2 da j

Отображение a : [— 1,1] ^ [0,1] играет ключевую роль в приведении интеграла (2) к виду (5): с одной стороны, оно выявляет структуру "плохих" функций задачи:

q = (1 — t)(1 + t)Q(t), и =1 — q = t2B(t), h- 1(a, s)Sda = £-1 IR-* 1(t)sin ^^L+A d t,

с другой же, переводит "подвижную" логарифмическую особенность из (2) в неподвижную - середину отрезка I в (5).

Для приближённой реализации интегралов в (5) используются квадратурные формулы, учитывающие специфику поведения подынтегральных функций. Их построение связано со способом аппроксимации функций (6) многочленами из подпространства Vn алгебраических многочленов степени не выше п — 1.

Пусть С[/] и С'[/] (к > 0 - целое) - соответственно пространство непрерывных и к раз непрерывно дифференцируемых на отрезке I = [—1,1] функций /; \\ f\\ = max \ f(t) \ -

чебышевская норма и En(f) = inf \\f — Нп\\ = \\f — Rn\\ (п > 0 - целое) - наилучшее

H„ev "

чебышевское приближение функции f Е С [I] многочленом R п из Vп.

Пусть ¿1, 12,... , tn - разбиение отрезка I несовпадающими узлами. Интерполяционный многочлен Pn(t; f) с этими узлами и функционал погрешности 5n(t, f) имеют вид:

п

Pn(t; f) = ^ f(ti)Um (t), 6n(t, f) = f(t) — Pn(t; f).

i=1

Здесь ип%(t) - фундаментальные многочлены лагранжевой интерполяции [2].

Какой точности приближения можно добиться, если заранее ограничить степень п приближающих полиномов и способ их построения? Чтобы найти ответ на этот вопрос, заметим, что для любого многочлена Hn(t) Е Vп справедливо тождество

Sn(t, /) = № - Pn(t; /) = № - Hn(t) - Pn(t; f-Hn). Из этого соотношения, применив неравенство Лебега [2], получаем оценку погрешности:

I Sn(t, /)| = | № - Pn(t; /)| < (1 + Лп) ВД),

где Лп - соответствующая константа Лебега. По теореме Лозинского-Харшиладзе [12] Лп > . На практике предпочтительнее разбиения отрезка с минимальным ростом константы Лп. Нули многочлена Чебышева первого рода Тп (t) = cos (п arccos t) (п > 1) вычисляются явно:

п (2 г - 1)

ti = cos -, U Е I, г = 1, 2,.. .п, п > 0.

2 п

Соответствующая им константа Лебега оценивается сверху по теореме Бернштейна [13]: Лп < 8 + ^ 1пп. При этом многочлен, интерполирующий функцию f(t) Е С[/], имеет вид

Рп(Ц /) = it f&) ТuTul)-f.), тп&) - (-1)*-1 п

<=i Tn(ti)(t - tiy V(1 - fc)(i + fc)'

а оценка погрешности удовлетворяет неравенству

I Sn(t, /)| = | f(t)-Pn(t;f)l< ^9 + 4 Inn) En(f). (9)

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

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

Хорошо известно [2,14], что числовой характеристикой En(f) обладает любая непрерывная функция f и, по теореме Вейерштрасса, lim En(f) = 0. Фундаментальное значение имеет существенное уточнение теоремы Вейерштрасса в форме так называемых неравенств Джексона: если f Е Ск[/] (к > 0), то

ж ак II/(к)||

En(/) < — min --—, а > 1 - абсолютная константа (10)

2 0<k<n пк

(усиленное1 неравенство Джексона [2, с. 307]).

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

Есть определённые преимущества в рассмотрении бесконечно дифференцируемых функций. Этот случай, возможно, и проще для исследования, но и ожидания здесь будут соответственно выше. Действительно, при f Е Сте[/], f Е ^n, II/II = ^(0) = 0,

1В статье [3] это неравенство приведено с опечаткой. Автор благодарен проф. Р.М. Тригубу, обратившему его внимание на это обстоятельство.

'(к)|| < G(k), lim ^йЩ = то определены следующие функции аргумента г Е [0, то):

к^-х

'G(0) при 0 < r< 1,

Л(0 = Г , G(k)

min

0<к<г г

к

при г > 1,

0

при 0 < < 1 ,

( ) = G( )

max < k | 1 < k < г и Л(г) = —^ > при г > 1,

причем выполнено

) . G(k) G[fl(r)] Л(г) = min —— = -—— и

0<к<г Гк

X'

(г)

— Rn|| =En(f) < 2 Л ( п) .

(11)

Теорема 1 [3]. При г > 1 функция в (г) целочисленна и неотрицательна, не убывает, непрерывна справа и стремится к бесконечности вместе с г; функция Л(г) строго монотонно убывает, непрерывна справа и стремится к нулю при г ^ то. Функция Л(г) терпит разрывы слева лишь в точках разрыва функции в (г). При этом для любого £ > 0 справедливо равенство

( )

Л( ) = Л( )

-Е и

Ti<X

> .

Здесь о0 = 0 и oi = ln Л(ri — 0) — ln Л(г\) для всех i > 0.

Следствие. Для любого р > 0 имеет место равенство limrг-Р Л(г) = 0.

Из неравенства (11) следует, что наличие информации о "большом" запасе гладкости функции f обретает, в силу теоремы 1, вполне конкретное и осязаемое на практике значение. С ростом п чебышевский аппроксимационный процесс самосовершенствуется, т.е. автоматически и непосредственно в дифференциальной природе черпает свою практическую эффективность, настраиваясь по фактической гладкости f на максимальный для данного п порядок сходимости в(п) (феномен ненасыщаемости [3]). Интуитивно такое свойство метода (9) означает, что с ростом "запаса" гладкости функции f скорость убывания погрешности метода к нулю возрастает. При этом потенциальные возможности развития метода зависят исключительно от скорости роста функции в(п) при увеличении параметра п. Максимума своей практической эффективности - экспоненциальной сходимости - метод обретает на классе Сх-гладких функций. Для f Е Сх[/] при условии

limk^x ^G(k)/ka < то, где а > 0, имеем | En(/)| < Се—@^п, здесь С, g - положительные постоянные.

Иначе говоря, конструкция метода (9) уже изначально такова, что способна вместить в себе, образно говоря, возможности бесконечного множества численных методов. Выявляя возможность "гибкого" участия запасов гладкости функции в оценке погрешности, теорема 1 устанавливает также, что приближенный метод, погрешность которого оценивается в терминах характеристик En(f), ненасыщаем [2,3]. В частности, таковым является метод

(9).

Свойство ненасыщаемости расширяет сферу практической применимости численного метода для класса Сх-гладких функций со следующим свойством.

Определение [3]. Функция f е Сх[I] обладает на отрезке I = [—1,1] пограничным слоем толщины е0, 0 < е0 < 1, если существуют такое малое положительное число т = т(е0) и такая положительная функция F(k), от е0 не зависящая, что при любом целом

к > 0 справедливо неравенство

Ifmw< I F(%lt) если I е n"[-1+т'1 -тЬ (12)

е0 kF(к), если t е I \ 1Т

где £0 = £0(т) > 0 - толщина пограничного слоя, а Р(к) - функция параметра к.

Специфика поведения функции £ Е С ^ [I] вблизи концов отрезка I отражается на сходимости к нулю с ростом п её чебышевских характеристик Еп(/). Наводящим соображением к этому служит:

Теорема (В.К. Дзядык [14]). Для того чтобы к-я производная функции / удовлетворяла на отрезке I = [— 1,1 ] условию Гёльдера с показателем 0 < а < 1, необходимо и достаточно, чтобы при любом целом п > к существовал такой алгебраический многочлен Qn Е Vп степени п — 1, что для всех Ь Е I справедливы оценки

Ак ( п-^ 1

| m -Qn(t)i< ^ (vr-*+£)

n k+" \ п /

где Ak — постоянная, не зависящая от t и п.

Приведённый результат усиливает неравенство Джексона (10), ибо при равномерной оценке En(f) < Мп~(k+a) устанавливает возможность приближения вблизи концов I c погрешностью 0(п ~2(k+a)). Придадим последнему наблюдению строгий смысл, связав его с определением пограничного слоя.

Теорема 2 [3]. Если функция f е С те[1] и выполнено условие (12), то

End) < 2 min ^^. (13)

2 0<k<n nk

Здесь коэффициенты Fk эффективно вычисляются по значениям F(к) из (12)(см. [3]).

Прикладное значение оценки (13) в том, что за счёт перераспределения пограничного слоя по всему отрезку I, толщину его удалось увеличить до значения у/~Ё0. При этом

гладкость функции f е Ck [I] характеризуется величиной £0 Fk, рост которой компенсируется, согласно теореме 1, выбором параметра п.

Действительно, при фиксированном п среди неравенств (13), отвечающим различным 0 < к < п, имеется наилучшее: его номер к0 = в(п) - это порядок той (максимальной) производной, которая участвует в формировании оценки (13). Производные более высоких порядков к > к0 могут, начиная с некоторого порогового значения птП = пт\п(ео), влиять на величину оценки (13) только в случае п > птщ. Таким образом, процедура нейтрализации пограничного слоя осуществляется, вследствие теоремы 1, за счёт выбора параметра п > пт{п: En(f) < ъ/2 mino<k<n £о-k/2Fk/пк < t, 0 < е < £о.

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

необязательно, чтобы f е Сте[!]: вполне можно обойтись конечным в(п) запасом произ-

_1/2

водных. Их количество при п ^ оценивается, как в(п) = 0(е0 / ).

4. Ненасьщаемые квадратурные формулы. Принимая во внимание, что

гс (¿) = ад г с) + ад Тс), ад = рп(ц гь) + ад гь)

и подставляя их в (5), получаем равенство

тл/ ч -1 ^ гс(U) [ Tn{t) - ^ TL(U) f Tn{t) / ^ч /1 л

rw=£ hш] ~tdt-e gtmj ]nltldt+iJn(s-r) (14)

где дп(в, Г) - функционал погрешности. Особые интегралы в (14) вычисляются явно, и на их основе в [15] построены квадратурные формулы:

+1 +1

/ f(t)dt = £ с,/ад pn ( A - /ад \t\dt = Y, и лад рП ( /). (15) i=1 1\ i=1

Здесь рП и рП - функционалы погрешности, а веса ск и Iк задаются формулами = 2 С(п, U) 2 L(n, U) (-1)г-1п

Ci = 2 rri,u \ , 1, = — 2 rri,U Ч , ТП W

ТП (ад г ТП (t„ ) ' пкг' ад-адад'

которые реализуются с помощью соотношений

[ ^ ]

{СП%}=±{ -1} " - -С- ] } (п> 0).

I Г ... Г 1 ~ (-1) J С[ ^ ]

1 ^ ]

При этом используются предварительно рассчитанные коэффициенты: С1 = 2, L1 = 2,

Ст = 1/(т - 0.5), Lm = -(т - 1.5)CmLm-1 - /(т - 1.5)/2, 2 <т <

п + 1

2

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

Теорема 3 [15]. Если п > 1 и нечетно, то веса с, и k в формулах (15) строго положительны, а их функционалы погрешностей рЩ, рП удовлетворяют неравенствам

\ рП(Л \< 4 Еп(f), \ рП( /) \ < 4 Еп(f). (16)

Из оценки (16) следует, что квадратурные формулы (15) ненасыщаемы; положительность же коэффициентов с,, I, обеспечивает им устойчивость к ошибкам округлений.

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

/ , ^ л \ к

< е(s) = 0.5ж\sin^s|, F(k) = е-к(s) max

f . 2 ж (a + s)d\

Г2 ^Т"Та)

(а, в) ЕКхК V 2 ¿а/

а за функцию /(а, в ) принять либо Гс (а, в), либо Гь (а, в) с учётом преобразования (4).

В результате у квадратурных формул (15) имеется реальная возможность для выполнения неравенств | рП (Гс)| < | рП (Гь)| < е, 0 < е < £0, которая и осуществляется,

согласно теореме 1, за счёт выбора параметра п > птщ(е0) в оценках (16). Это влечёт, исходя из требуемого количества к0 = в(п) производных у функций Гс(1), Гд(1) Е С2р+1[1], нейтрализацию пограничного слоя в (5). Выбор параметра р подчиним условию нейтрализации пограничного слоя, т.е. неравенству р > [-^— ].

В прецизионных алгоритмах численного решения осесимметричных краевых задач следует использовать "гибкие" методы вычисления полных эллиптических интегралов (например, указанные в [10]), поскольку именно они адекватно соответствуют природе рассматриваемых вопросов и удачно завершают конструкцию ненасыщаемого алгоритма численной аппроксимации интеграла Гаусса (1).

5. Алгоритм вычисления. Численные примеры

Компьютерная реализации интеграла Гаусса (1) осуществляется по следующим формулам:

(п п \

^ С i Г с Цг) + гГь(1г)\ + Qn( s, Г), е = 0.5п sinns. (17)

г=1 г=1 J

При этом функции Гс (t), Ti,(t) заданы формулами (6), а веса с i, I ¿ - формулами (15). Величины a i = а§(tí, s) (1 < г < п) восстанавливаются по фиксированному значению s Е (0,1) и заданному упорядоченному массиву чисел {t¿ = cos ж(22П-1}:

. (Ji-S) . (ai + S ) . t i = Sin---/ Sin---, 1 < г < П.

Восстановление a i производится методом Ньютона посредством следующего алгоритма:

k+i к к к

У =У - f (у)/f '(у), к = 0, 1, 2,... - номер итерации,

tf \ 2 ■ U • ^(х + s) \ , ticoS

j(x) = х - s--arcsm (ti sin-1, j (х) = 1--. =.

^ -t2 sin2

0 о

Вопрос о выборе начального приближения У=ai (1 < % < п) в методе Ньютона важен. В рассматриваемом случае он решается достаточно просто:

о , -i Л + 4. n s\ о n(t i+i -t i) . 2n(ai + s) i

Ji =1 -n 1 1 - ti ctg — , a i+i = л + 2-;- sin2---, 1 < г <п - 1.

V 2 J n sinns 2

Указанный итерационный процесс обладает квадратичной сходимостью: уже 2-3 итерации обеспечивают не менее 10 точных десятичных разрядов для s Е (0.01, 0.99 ). Нетрудно

указать быстрый итерационный процесс и для s Е (0.001,0.999): при этом следует лишь

о о

несколько модифицировать выбор значений = a (1 < < п).

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

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

Ограничимся рассмотрением эллипсоидов вращения с меридиональным сечением:

7( 5 ) = {а вткз, ЬооБкв}, к = Ь/а, вЕ [0,1],

отдав предпочтение эллипсоидам, вытянутым вдоль оси симметрии г (Ь > а), как наиболее "трудным" с точки зрения проведения численных расчетов областям.

Интеграл Г( в) вычислялся в точках Sj = ]/100 (1 < ] < 100) по формулам (17).

В табл. 1-5 указана разность к — Г(Sj) - абсолютная ошибка - в первых десяти точках вблизи полюса 7о: именно в них влияние пограничного слоя наиболее ощутимо, хотя его присутствие, конечно, проявляется в каждой расчётной точке Sj (1 < ] < 100), и чем больше к, тем оно ощутимее. При этом рассмотренный диапазон удлинений 1 < к < 100 демонстрирует реальные возможности ненасыщаемого численного метода.

Смысл результатов, представленных табл. 1-5, состоит в следующем: если подынтегральные функции в (5) при р > р0 существенно более гладкие, чем при р = р0, то при прочих равных условиях метод вычисляет Г(вj) с точностью большей, чем в случае р) = р0. Случай р) = 4 соответствует известному [16] способу аппроксимации полных эллиптических интегралов. Добавим к этому, что при расчёте табл. 5 одна только замена в алгоритме параметра р) = 5 нар =10 увеличивает точность до 8 десятичных разрядов. Другие примеры использования построенной методики можно найти, например, в работе [11].

Автор благодарен проф. В.Л. Васкевичу за помощь в редактировании статьи.

Таблица 1. Интеграл Гаусса Г(б)

а = Ь = 1/к, Sj = j /100, п = 40

р 4 10 20 40

3 к — Г( Sj ) к — Г( Sj ) к — Г( Sj ) к — Г( Sj )

1 0.00074996015 0.00085952336 0. 00094117975 0. 00096930719

2 0.00003620611 0.00005709499 0. 00005785123 0. 00005327380

3 0.00002540592 0.00000476493 0. 00000401157 0. 00000431754

4 0.00002168817 0.00000057541 0. 00000041157 0. 00000041578

5 0.00002993194 0.00000009031 0. 00000021862 0. 00000023572

6 0.00003338771 0.00000014374 0. 00000000286 0. 00000001958

7 0.00003562241 0.00000017582 0. 00000001728 0. 00000000286

8 0.00003693635 0.00000016574 0. 00000000166 0. 00000001053

9 0.00003752035 0.00000016399 0. 00000000158 0. 00000001187

10 0.00003751165 0.00000016414 0. 00000000009 0. 00000000887

Таблица 2. Интеграл Гаусса Г (б)

а = Ь = 1/тг, sj = j /100, п = 80

р 4 10 40 80

3 т - г(в]) т - г(^) т - г(^) т - г(^)

1 0.00000791729 0.00000030422 0.00000036784 00000034329

2 0.00001499472 0.00000008073 0.00000000447 00000000455

3 0.00002152174 0.00000011099 0.00000000197 00000000200

4 0.00002704800 0.00000013839 0.00000000216 00000000217

5 0.00003145025 0.00000015930 0.00000000203 00000000203

6 0.00003471950 0.00000017262 0.00000000045 00000000045

7 0.00003693750 0.00000018305 0.00000000148 00000000148

8 0.00003823076 0.00000018725 0.00000000105 00000000107

9 0.00003875586 0.00000018790 0.00000000094 00000000093

10 0.00003866926 0.00000018571 0.00000000101 00000000080

Таблица 3. Интеграл Гаусса Г (б)

а = Ь = 1/тг, Sj = j /100, п = 100

Р 4 10 50 100

3 т - г(в]) т - г(в]) т - г(в]) т - г(в])

1 0.00000775635 0.00000005011 00000000742 00000000831

2 0.00001508046 0.00000008207 00000000443 00000000443

3 0.00002164322 0.00000011590 00000000492 00000000490

4 0.00002718896 0.00000013937 00000000080 00000000079

5 0.00003160697 0.00000016184 00000000192 00000000190

6 0.00003488173 0.00000017458 00000000041 00000000041

7 0.00003710023 0.00000018590 00000000148 00000000147

8 0.00003838814 0.00000018960 00000000059 00000000061

9 0.00003890500 0.00000018984 00000000018 00000000015

10 0.00003880950 0.00000018815 00000000080 00000000080

Таблица 4. Интеграл Гаусса Г^)

а = Ь = 1 /к, Sj = j /100, п = 100

р 5 25 55 95

3 к — Г( Sj ) к — Г( Sj ) к — Г( Sj )

1 0.00000309766 0.00000000860 0.00000000742 0.00000000827

2 0.00000600681 0.00000000445 0.00000000443 0.00000000443

3 0.00000861302 0.00000000489 0.00000000491 0.00000000490

4 0.00001080007 0.00000000080 0.00000000082 0.00000000080

5 0.00001253228 0.00000000191 0.00000000192 0.00000000192

6 0.00001379995 0.00000000039 0.00000000040 0.00000000042

7 0.00001464497 0.00000000144 0.00000000148 0.00000000147

8 0.00001511599 0.00000000061 0.00000000061 0.00000000061

9 0.00001528089 0.00000000018 0.00000000017 0.00000000016

10 0.00001520526 0.00000000080 0.00000000081 0.00000000080

a = 1, b =

Таблица 5. Интеграл Гаусса r(s)

100, Sj = j /100, n = 2000, p = 5

3 п - Г(Sj) 3 п - Г(Sj)

1 0.00000019329 26 0.00000269393

2 0.00000123943 27 0.00000225247

3 0.00000186769 28 0.00000165033

4 0.00000199312 29 0.00000201233

5 0.00000196848 30 0.00000241428

6 0.00000207644 31 0.00000285311

7 0.00000201963 32 0.00000218227

8 0.00000238756 33 0.00000135977

9 0.00000224457 34 0.00000157173

10 0.00000206489 35 0.00000179447

11 0.00000189213 36 0.00000202522

12 0.00000234309 37 0.00000226069

13 0.00000216461 38 0.00000249735

14 0.00000224396 39 0.00000273245

15 0.00000215630 40 0.00000296148

16 0.00000191624 41 0.00000318214

17 0.00000214350 42 0.00000338972

18 0.00000234187 43 0.00000358196

19 0.00000171960 44 0.00000111728

20 0.00000243399 45 0.00000116935

21 0.00000236408 46 0.00000121328

22 0.00000217242 47 0.00000124865

23 0.00000287342 48 0.00000127462

24 0.00000165453 49 0.00000129018

25 0.00000213453 50 0.00000129555

СПИСОК ЛИТЕРАТУРЫ

1. V.N. Belykh To the problem of evolutionary "blow-up" of axially symmetric gas bubble in ideal incompressible fluid (main constructive hypothesis) // Proceedings of Intern. Conf. dedicated to M.A. Lavrentyev on the occasion on his birthday centenary. Ukraine, Kiev, 2000, p. 6-8.

2. Бабенко К.И. Основы численного анализа. М.: Наука, 1986, 744 с. (2-е издание: М.; Ижевск: РХД, 2002. 847 с.).

3. Белых В.Н. О свойствах наилучших приближений С^-гладких функций на отрезке вещественной оси (к феномену ненасыщаемости численных методов) // СМЖ. 2005. Т. 46, № 3. С. 483-499.

4. Агмон С., Дуглис А., Ниренберг Л. Оценки решений эллиптических уравнений вблизи границы. ИЛ., М., 1962. 205с.

5. Белых В.Н. Внешняя осесимметричная задача Неймана для уравнения Лапласа: ненасыщае-мые методы численного решения // ДАН. 2007. Т. 417, № 4. С. 442-445.

6. Н.Н. Анучина, К.И. Бабенко, С.К. Годунов и др. Теоретические основы и конструирование численных алгоритмов задач математической физики. М.: Наука, 1977.

7. Васкевич В.Л. Гарантированная точность вычисления многомерных интегралов // Дис. ... д-ра. физ.-мат. наук. Ин-т математики им. С.Л. Соболева СО РАН. Новосибирск, 2003. 243 с. РГБ ОД, 71:04-1/275.

8. Рамазанов М.Д. Новый алгоритм асимптотически оптимальных решётчатых кубатурных формул // Уфимский математический журнал. 2010. Т. 2, № 3. С. 63-82.

9. Гюнтер Н.М. Теория потенциала и ее применение к основным задачам математической физики. М.: Гостехтеоретиздат, 1953. 415c.

10. Белых В.Н. Алгоритмы вычисления полных эллиптических интегралов и некоторых свя-зянных с ними функций // Сиб. журнал индустр. матем. 2012. Т. 15, № 2 (50). С. 21-32.

11. Белых В.Н. К проблеме обтекания осесимметричных тел большого удлинения потоком идеальной несжимаемой жидкости // ПМТФ. 2006. Т. 47. № 5. С. 56-67.

12. Даугавет И.К. Введение в теорию приближения функций. Л., Изд-во Ленингр. ун-та, 1977. 184 с.

13. Бернштейн С.Н. Собрание сочинений. Т.2, АН СССР, М., 1954. 627с.

14. Дзядык В.К. Введение в теорию равномерного приближения функций полиномами. М.: Наука, 1977.

15. Белых В.Н. Алгоритмы без насыщения в задаче численного интегрирования // ДАН. 1989. Т. 304, № 3. С. 529-533.

16. Абрамовиц М., Стиган И. Справочник по специальным функциям. М.: Наука, 1979. 832 с.

Владимир Никитич Белых,

Институт математики им. С.Л. Соболева,

Проспект акад. Коптюга, 4,

630090 г. Новосибирск, Россия

E-mail: belykh@math.nsc.ru