О вычислительной устойчивости одного метода асимптотической стабилизации Текст научной статьи по специальности «Математика»

9. Gathen J. von zur, Gerhard J. Modern computer algebra. Cambridge: Cambridge University Press, 1999.

10. Серпинский В. 250 задач по элементарной теории чисел. М.: Просвещение, 1968.

11. Rosser J., Schoenfeld L. Approximate formulas for some functions of prime numbers // 1ll. J. Math. 1962. 6. 64-94.

Поступила в редакцию 06.03.2006

УДК 519.6

О ВЫЧИСЛИТЕЛЬНОЙ УСТОЙЧИВОСТИ ОДНОГО МЕТОДА АСИМПТОТИЧЕСКОЙ СТАБИЛИЗАЦИИ

А. А. Корнев, А. В. Озерицкий

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

Пусть в банаховом пространстве Н с нормой \\ ■ \\ задано достаточно гладкое отображение 5 : Н — Н. Обозначим через № -(3, О) устойчивое многообразие подмножества Ос Н:

№- (3, О) = {то еО: Зт+г е О, ш^г = 3(шг), г = 0,1,2,...} .

Отметим, что устойчивое многообразие состоит из множества точек то, всегда остающихся в О при действии оператора 5. Если О является окрестностью неподвижной (периодической) точки либо окрестностью траектории гиперболического типа, то структуру устойчивого многообразия можно уточнить. Соответствующее утверждение принято называть обобщенной теоремой Адамара-Перрона (см. работу [1] и имеющийся в ней обзор).

Рассмотрим задачу проектирования заданной точки ао С О на многообразие №- в окрестности О нулевой неподвижной точки 5(0) = 0.

Будем считать, что для оператора 5 существуют операторы проектирования Р+,Р- : Н — Н, ограниченный линейный оператор Ь, непрерывное отображение Я(и) = 3(и) — Ьи, такие, что имеют место следующие (см. [1]) условия гиперболичности (А):

Р+ + Р- = I, \\Р+1|, \\Р-1| < С; Ь(Р+Н) = Р+Н, Ь(Р-Н) С Р-Н; \\Ьь\\ > ^+\М\, Уь е Р+ Н, > 1; \\bw\l < ^^Ц, Уw е Р-Н, ^- < 1; ЦЩиг) — Е(п2)\\ < 0(т&х{\\и1\\, \\п2\\})\\и1 — и2\\, Ущ е Н

с непрерывной положительной неубывающей функцией 9(-), 9(0) = 0. Так как подпространство Р+Н обычно конечномерно, а Р- Н имеет конечную коразмерность, то операторы проектирования Р± удобно задавать следующими наборами базисных векторов: Р+Н =< в+,. ..,в+ >, Р^Н =< е— ,...,е- >, \\\ = 1 для г = 1,...,%о.

Изложим основанный на работах [1-3] метод численного построения многообразия №Запишем оператор 3(и) = Ьи + Я(и) для и = V + w, V е Р+О, w е Р-О в виде 3(и) = 3+(и) + 3-(и), где 3±(и) = Р±3(и). При этом

3+(ь + w) = Ь+ь + Я+(ь + w), Ь± = Р±Ь, 3-(ь + w) = Ь^ + Я-(ь + w), Я±(-) = Р±Я(-).

Будем искать [2] устойчивое многообразие в виде

№-(О) = {Р-и + /(Р-и), и еО}

с некоторой функцией f (w) : P-O ^ P+ O. Выпишем условие инвариантности устойчивого многообразия W- относительно оператора S:

P+Sf (w) + w) = f^P-S (f (w)+w)),

или в эквивалентном виде:

P+L(w + f (w)) + P+ R(f (w) + w) = f(P-L(w + f (w)) + P-R(f (w) + w)).

Из условий (А) следует инвариантность подпространств P±H относительно оператора L и обратимость L на подпространстве P+H. Поэтому данное соотношение можно переписать (см. [3]) следующим образом:

f (w) = L-1(f (P-S(w + f (w))) - P+ S(w + f (w))) + f (w), (1)

а для решения полученного уравнения f (w) = F(f (w)) относительно неизвестной функции f применить итерационный процесс

fk+i(w) = F(fk(w)), fo = 0. (2)

В работах [1, 2] доказана теорема о сходимости в некоторой окрестности O итерационного процесса (2) к устойчивому многообразию W-. Рассмотрим вопрос об устойчивости данного алгоритма относительно ошибок вычисления операторов P±. Во-первых, отметим, что вывод соотношения (1) существенно опирается на точную инвариантность подпространств P±H относительно оператора L. Во-вторых, условия (А) гарантируют обратимость и сжимаемость оператора L только на подпространстве P+H, что также принципиально при построении доказательства. На всем пространстве H оператор L может быть плохо обусловлен либо вырожден, поэтому формально даже малые ошибки в построении P±H недопустимы. При практической реализации итерационного процесса (2) операторы проектирования P± вычисляются с некоторой погрешностью (например, с машинной точностью), влияние которой на окончательный результат f(w) может оказаться катастрофическим.

Сформулируем метод проектирования на многообразие W- в случае приближенно вычисленных векторов e±, i = 0,...,io. Пусть оператор L невырожден и известны такие операторы проектирования P±, что выполнены следующие условия (А):

\\Lu\\ < C+\\u\\, \\L-lu\\ < C-\\u\\, Уп e H; P+H =< ё+,..., e+ >, P-rH =<ё-,..., e- >,

e± = e± + \\e±\\ < Ф±\ i = 1,...,io.

Рассмотрим класс В7 (O) всех непрерывных отображений f (w) : P-O ^ P+O, таких, что

f(0) = 0, \\f(wi) - f(w2 )\\ < Y\\wi - w2 \\, Y = const < то.

В пространстве B7 (O) зададим норму \f \ = max \\f (w)\\.

weP-O

Имеет место следующая теорема.

Теорема. Пусть выполнены условия (А), (А). Тогда найдутся такие £,Y,r > 0, что в окрестности O = {u : \\P-u\\ < r, \\P+u\\ < Yr] итерационный процесс

fk+i(w) = P+L~l(f k(P-S(w + fk(w))) - P+S(w + fr(w))) + f k(w) (

сходится в пространстве B7 (O) со скоростью геометрической прогрессии с любого начального приближения fo e В^ (O). Предельная функция является решением уравнения (1). Для каждого элемента w eO выполняется условие \\Sl(w + f (w))\\ dO, i = 0,1,..., и верна оценка

\\S*(w + ёН)\\ < C ?\\w + ёН\\, C = const < ос,ё< 1.

вестн. моск. ун-та. сер. 1, математика. механика. 2007. № 1

35

Замечание. В работах [3, 4] метод (3) применялся с целью повышения вычислительной устойчивости алгоритма (2).

Данная теорема может быть доказана на основе метода сжимающих отображений по схеме, предложенной в работе [2]. Наличие векторов ошибки ё± приводит к появлению в расчетных формулах дополнительных слагаемых, зависящих от е, С±. Однако при условии малости величины е можно показать, что на каждом шаге алгоритма (3) функция f принадлежит классу В7 (О), а соответствующий оператор перехода ¡к+1 = Г[/к) является сжимающим в указанной норме. Это позволяет доказать сходимость итерационного процесса. Завершает доказательство проверка, что каждая неподвижная точка отображения (3) является решением уравнения (2).

Рассмотренные методы (2), (3) построения отображения f являются основой для решения задачи проектирования на устойчивое многообразие вдоль подпространства допустимых смещений С =< 1г, ...,1^0 >, а также в окрестности траектории (см. [3, 4]). Соответствующие алгоритмы также могут применяться для стабилизации по начальным данным неустойчивых решений нестационарных уравнений в частных производных.

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

^ + .1(ф, Аф) + .1(ф, 1 + к)+ сгАф - А2ф = 0, где ф(Ь,ф,Л) — функция тока, ф е [0,п/2] — широта, Л е [0, 2п] — долгота; А — оператор Лапласа на

единичной сфере; J(u, v) =

cos ф

ди ду _ ди dv_

д\ дф дф д\

l = 20 sin ф — параметр Кориолиса, 0 = 7, 27 ■ 10

-5.

h = 1, 75 cos 2Х sin2 2ф — ортографические неоднородности подстилающей поверхности; а — коэффициент приземного трения. Если значение а > 0, то данное уравнение (см. [5]) описывает динамику атмосферы над двумя континентальными массами, разделенными двумя океанами; высоты континентов понижаются к полюсам и к экватору. Возьмем а = -20. Формально выбор а < 0 не является физически правильным, однако на малых временных отрезках позволяет качественно моделировать сильную неустойчивость траекторий данного уравнения.

Рассмотрим решение задачи для случая полусферы (северного полушария) c начальным условием ■0|t_o = v и краевыми условиями непротекания на экваторе §х1(/>-о = о =

При построении конечно-разностной аппроксимации применим схему Аракавы по пространству и Кранка-Николсона по времени. Приведем результаты расчетов проектирования на устойчивое многообразие для следующей функции v = uo, которая является заведомо неустойчивым начальным условием данной задачи:

Uo = ¡~2 при| <Л<§, f <0<f;

Icos Л sin 2ф иначе.

Построим оператор L аналитически по разрешающему оператору S соответствующей разностной задачи. Подпространства P± H вычислим методом Арнольди как инвариантные подпространства оператора L. При этом подпространство P+H соответствует собственным значениям, по модулю большим единицы, а подпространство P- H — меньшим единицы.

По заданной функции uo найдем такую функцию l gCc. H, что uo +1 G W-(S, O). Базис li подпространства допустимых смещений L построим по вычисленным базисным функциям e+ подпространства P+H следующим образом:

li(фn, Лт) —

о, Фе[f,f], Ле [0,2тг];

e+ (ф n, Лт) иначе.

Отметим, что такой выбор базиса означает, что функция ио не изменится в прилегающей к экватору полосе. На рис. 1, 2 изображены соответственно начальная функция ио и найденная функция I при следующих значениях сеточных параметров: 32 х 24 (долгота х широта); шаг по времени 0, 01.

В таблице приведены результаты численных расчетов с начальными условиями V = ио, ио + I до момента Т, где Т — время стабилизации решения, т.е. время, в течение которого норма решения с соответствующим начальным условием V убывает.

V Т №,т)|| IMI

ио 0,01 16,75 2,0

ио +1 0,67 0,0085 2,23

1

Рис. 1 Рис. 2

В общем случае аналитическое выделение оператора L является технически сложной задачей. Поэтому представляет интерес возможность приближенного построения линейной части L на основе численных расчетов с оператором S и последующего нахождения P± по оператору L. Это соответствует приближенному вычислению операторов проектирования для оператора L. В результате такого построения подпространств P± и дальнейшего решения задачи проектирования по методу (3) была найдена функция I ЕС, практически не отличающаяся от исходной: \\l — l Ус = 1, 7 ■ 10-6. Имеющиеся различия обусловлены точностью решения задачи (3) и сильной неустойчивостью разрешающего оператора S данной задачи.

Таким образом, существенные ошибки в вычислении проекторов практически не сказались на точности расчетов, однако повлияли на скорости сходимости алгоритма. Для данной задачи количество итераций, необходимое для нахождения значения функции f (w) с точностью 10-12 в случае проектирования с помощью операторов P±, равнялось 11, а с помощью операторов P± равнялось 22.

Отметим, что полученная теорема показывает устойчивость итерационного процесса (3) для решения задачи проектирования на многообразие W-, хотя оператор S исходных уравнений является существенно неустойчивым.

Авторы признательны А. С. Грицуну и Ю. М. Нечепуренко за полезные обсуждения.

Работа выполнена при поддержке гранта РФФИ № 05-01-00864 и гранта NWO 047.016.008.

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

1. Аносов Д.В. Геодезические потоки на замкнутых римановых многообразиях отрицательной кривизны // Тр. Матем. ин-та АН СССР. 1967. 90.

2. Ладыженская О.А., Солонников В.А. О принципе линеаризации и инвариантных многообразиях для задач магнитной гидродинамики // Зап. науч. сем. ЛОМИ. 1973. 38. 46-93.

3. Корнев А.А. Об итерационном методе построения "усов Адамара" // Журн. вычисл. матем. и матем. физ. 2004. 44, № 8. 1346-1355.

4. Корнев А.А. Метод асимптотической стабилизации по начальным данным к заданной траектории // Журн. вычисл. матем. и матем. физ. 2006. 46, № 1. 37-51.

5. Дымников В.П., Филатов А.Н. Основы математической теории климата. М.: ВИНИТИ, 1994.

Поступила в редакцию 28.04.2006