Градиентные методы решения задач Стокса Текст научной статьи по специальности «Математика»

ISSN 2074-1871 Уфимский математический журнал. Том 8. № 2 (2016). С. 22-38.

УДК 517.9

ГРАДИЕНТНЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ СТОКСА

И.И. ГОЛИЧЕВ, Т.Р. ШАРИПОВ, Н.И. ЛУЧНИКОВА

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

Ключевые слова: задача Стокса, оптимальное управление, градиентный метод, конечно-разностная схема.

Mathematics Subjects Classifications: 49M20, 35Q30, 93C05

1. Введение. В ограниченной области П С Rn (п =2, 3) с гладкой границей S G С2 (т.е. дважды непрерывно дифференцируемой) рассматривается задача Стокса

- vAv = f -Vp, v\s = 0, (1)

div v = 0, (2)

(P, iWn) = ° (3)

где f = (f{1),...,f(n)) G ¿2(П) = (L(П)Г, v = (w(1),...,v(n)) - вектор скорости, p -давление, v > 0 - коэффициент кинематической вязкости.

Согласно [1, гл. III, §5, теорема 2], в данных условиях задача (1)-(3) имеет единственное сильное решение v G H2(П) = (Н2(П))п, Vp G Х2(П). Здесь и далее используются стандартные обозначения для пространств Соболева Нг(П) = W2 (П), I = 1, 2 ...

Будем рассматривать задачу (1) - (3) как обратную задачу к задаче (1), (3), в которой неизвестен градиент давления Vp, но задано дополнительное условие (2). Такую задачу можно сформулировать как эквивалентную исходной задачи оптимального управления:

J(u) = 1 lldivv(u)|||2(Q) ^ inf, u G G(fi),

где G(Q) = {u G L2(n) : u = Vp, p G H 1(П), (p, 1)l2(q) = 0} - градиентная составляющая в ортогональном разложении L2(H), а v(u) - решение системы

- vAv(u) = f- u, v(u)|s = 0. (4)

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

I.I. Golichev, T.R. Sharipov, N.I. Luchnikova, Gradient methods for solving Stokes problem. © ГоличЕв И.И., Шарипов Т.Р., ЛучниковА Н.И. 2016.

Поступила 9 декабря 2015 г.

Наряду с задачей минимизации функционала J(u) в пространстве L2(H) также будем рассматривать задачи минимизации J(Vu) в пространствах Н^П) и ¿2(П), в которых за управление и берется само давление р. Поскольку градиенты функционала в этих случаях имеют различный вид, нам удобно их обозначить в соответствии с рассматриваемыми в работе задачами:

Задача I. Найти минимум функционала J0(u) = ||divv(u)|||2(q)/2 на множестве U0 = G(H), где v(u) - решение задачи (4).

Задача II. Найти минимум функционала J1(u) = ||divv(w)|||2(п)/2 на множестве U1 = {и Е Н^П) : (и, 1)l2(q) = 0}, где v(u) - решение задачи

- и Av(u) = f -Vu, v(u)|s = 0. (5)

Задача III. Найти минимум функционала J2(u) = ||divv(w)|||2(П)/2 на множестве U2 = [и Е L2(Q) : (и, 1)l2(q) = 0}, где v(u) - обобщенное решение задачи (5).

Решение задач I-III будем искать методом проекции градиента [6], [7]:

uk+i = Pui (ик - акJ'(uk)), I = 0, 1, 2, (6)

где Риг - оператор проектирования на множество Ui, а J/ (ик) - градиент функционала Ji (ик) в точке ик 1.

В дальнейшем изложении будет показано, что

J0 (u) = w(u), (7)

Ji (и) = p(w(u)), (8)

J2 (u) = - div w(u), (9)

где w(-) - сопряженное состояние системы (4), определяемое для всех задач как решение задачи

- ^Aw(-) = Vdivv(^), w(^)|s = 0, (10)

а p(w) определяется из ортогонального разложения вектора w(u) = Vp(w) + (р на градиентную и соленоидальную составляющие.

Для всех рассматриваемых задач найдем операторы проектирования и покажем, что задачи I и II эквивалентны.

Отметим, что в [4] рассмотрен аналог задачи II для обобщенной задачи Стокса, а для построения итерационного процесса использована общая теория сопряженных операторных уравнений. Тестовый расчет из этой работы используется нами для верификации расчета по формуле (6). К близким по структуре дифференциальных итерационных процессов настоящей работы приводит также метод неполного и полного расщепления граничных условий [5], где применяется теория граничных операторов Пуанкаре-Стеклова.

Далее везде через v(u) обозначается решение задачи (4) при f = 0, а через v(и) - решение задачи (5) при f =0.

2. Дифференцируемость функционала J0 (u). Обозначим через L оператор на множестве Н"^(П) = jz Е Н2(П) : z|s = 0} равенством L(z) = -vAz. Заметим, что область значений оператора L совпадает со всем пространством L2(П). Действительно, пусть f - произвольный элемент из L2(H); тогда задача

L(z) = f (11)

имеет, согласно [8, гл. II §7, теорема 7.1], единственное решение из Н2(П).

Используя второе энергетическое неравенство [8, гл. II, §6, формула (6.29)]2

||z||H2(n) < Ci||Az||£2(fi), (12)

1В целях обобщенной записи мы не выделяем жирными шрифтом вектор и для задачи I.

2Если область И выпуклая, то константа С1 равна единице.

ние

с учетом очевидного для задачи (11) равенства ||Az||^2(п) = V ЧШНх^п), получаем используемую далее оценку

Н^НЧп) < С1и-1Щ12(П). (13)

Пусть и и Ь - произвольные элементы из Щ. Учитывая определение оператора Ь, нетрудно видеть, что у(и + Ь) — у(и) = V(Ь), где V(Ь) = Ь-1 (—Ь). Тем самым,

Л (и + Ь) — Л (и) = 2 ||&У У(и + Ь)У|2(П) — 2 ||&у у(и) || ¿2 (п) =

= (ё1уу(и), а1у V(Ь))Ь2(П) + 1 ||ё1у V(Ь)||2(П). (14)

Преобразуем первое слагаемое в правой части последнего равенства с использованием интегрирования по частям, условия V(Ь)|^ = 0 и самосопряженности оператора Ь:

(&У v(u), V(Ь))^2(п) = — (У^У v(u),L-1(—Ь))£2(п) = v(u), Ь)^(п).

Таким образом, с учетом обозначения '(и) = L-1Vdiv v(u) тождество (14) принимает вид

Ми + Ь) — 7о(и) = ('(и), Ь)£2(п) + 2Иу V(Ь)||2(п). (15)

Умножая левую и правую части (11) на z и интегрируя по частям, получаем соотноше-

12

откуда, и из неравенства Фридрихса

|М|£2(П) < со^Н^п), (16)

имеем:

|^||£2(П) < со^Ш^п). (17)

Из (17) и легко проверяемого неравенства

Иуz|L2(п) < ^^^(п) (18)

следует, что

||divV(Ь)|Ь2(п) < ^псо V-1 ||Ь||£2(п), (19)

поэтому главная линейная часть приращения функционала </о(и) определяется выражением ('(и), Ь)£2(п).

Покажем, что градиент ^(и) удовлетворяет условию Липшица. Пусть и(1), и(2) е ио, а '(1), '(2) - соответствующие им решения задачи (10). Заметим, что

Д'(1) — '(2)) = Vdiv v(u(1)) — Vdiv v(u(2)) = Vdiv V (и(1) — и(2)), поэтому из неравенств (16), (17) получаем

||'(1) — '(2)|£2(П) < с2 ^-1||VdivV(и(1) — и(2))|£2(п). (20)

Учитывая то, что V(и(1) — и(2)) является решением уравнения (11) с правой частью Ш = —(и(1) — и(2)), для которого справедлива оценка (13), из (20) получаем:

К (и(1)) — (и(2))||£2(П) = ||'(1) — '(2) ^(П) < ^ ^ (и(1) — и(2))||н 2(п) <

С1 ^-2||и(1) — и(2)|£2(п). (21)

Таким образом, доказана

Теорема 1. Пусть f е ^(П), где П - ограниченная область с границей Б е С2; тогда функционал (и) дифференцируем по Фреше на и0, его градиент определяется по формуле (7) и удовлетворяет условию Липшица с константой Ь0 = у/пс° С1 и-2, где со и с1 - константы из неравенств (16), (12).

3. Дифференцируемость функционала ,]1 (и) На множестве и1 введем метрику, эквивалентную метрике пространства Н ^П), по скалярному произведению (М)я01(п) = (Va, УЬ)£2(П).

Пусть и и к - произвольные элементы из 171. Аналогично (15) легко видеть, что

Зг(и + к) - Ш = (ы(и), Ук)^(п) + 1 Иу V(к)||2(п), V(к) = £-1(-Ук), (22)

причем, с учетом разложения вектора ы(и) = Рс(п)ы(и) + ф на градиентную и соленои-дальную (diy ф =0) составляющие,

(ы(и), Ук)^(п) = (РС(п)™(и), Ук)^(п). (23)

Далее известно [9], что Vd € £2(П) = Ур, где р есть решение задачи Неймана

Ар = div d,

оп

= (d ■ п)|5. (24)

в

Из (22), (23), (24) следует:

31(и + к) - = (рЫ,к)Я0(п) + 1 ||diy ^ (к)||2(П)

где р(ы) - решение задачи

—р(ы)

Ар(ы) = div ы(и),

—п

= (ы(и) ■ п)|5. (25)

я

Аналогично (19) имеем оценку

Иу V(к)|к(п) < Со ^-1|ук|£2(п) = у/пСо т^-11 к|я0(п),

откуда следует, что главная линейная часть приращения функционала 3\_ (и) определяется выражением (р(ы), к)Я0(п).

Покажем, что градиент (и) удовлетворяет условию Липшица. Пусть м(1) и м(2) принадлежат [Д, а ы(1) и ы(2) - соответствующие им решения задачи (10); тогда

К (и(1)) - 3[ (И(2))|Я1(П) = ||р(«(1)) - р(«(2))||я1(П) = ||Ур (и(1)) - Ур(и(2) )|£2(П) =

= ||^Ъ(п)Ы(1) - РС(п)Ы(1)|£2(п) < ||Ы(1) - Ы(2)||£2(п).

Учитывая первое из неравенств (21), а также то, что V(м(1) - м(2)) является решением уравнения (11) с правой частью f = -У(и(1) - и(2)) и принимая во внимание оценку (13), получим

|№(1)) - ^(«(2))||я01(П) < ||Ы(1) - Ы(2)||£2(П) < у/п с2с1 ^-2|У^(1) - У«(2)|Ь2(П) =

= Vп с0с1 и-2Ци(1 -и(2)ЦН0(п).

Теорема 2. Пусть выполнены условия теоремы 1; тогда функционал /1(м) дифференцируем по Фреше на и1, его градиент определяется по формуле (8) и удовлетворяет, условию Липшица с константой Ь0.

Если в задачах I и II начальные приближения в итерационных процессах (6) при I = 0 и 1 = 1 связаны равенством и0 = Уи0, то щ = Уик при всех к = 1, 2,...

Доказательство. Действительно, если щ = Уик, то с очевидностью v(ufc) = v(uk) и ) = ы(ик), поэтому

ик+1 = Р0(п) (ик - ак30(ик)) = ик - акРС(п)ы(ик) = = ик - акУрк(ик) = У(ик - акрк(ик)) = Уик+1. Таким образом, задачи I и II можно считать эквивалентными. □

4. Дифференцируемость функционала (и). Пусть и и к - произвольные элементы из и2. Для доказательства формулы (9) воспользоваться непосредственно, как в предыдущих разделах, интегрированием по частям здесь невозможно, поскольку не гарантирована принадлежность функций v(u) и '(и) пространству Н2(П). Воспользуемся предельным переходом; выберем последовательности {ип}, {кп}, содержащиеся в и1 так, что ип ^ и, кп ^ к в Ь2(П) при п ^ го.

На последовательностях {ип}, {кп} справедливо равенство (22):

■Ыип + кп) — Мип) = ('(ип), Vкп)i2(п) + 1 ||divV(ЮЦ^п). (26)

Покажем, что div'(ип) ^ div'(и), divV(кп) ^ divV(к) при п ^ го в метрике Ь2(П).

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

V Н^Ьф) = № ^2(П) + (Щ div

применение неравенств (18) и (16), к которому дает неравенство

|^||£2(П) < ^ЫШ^п) + ^\\и\\ЫП)). (27)

Разность v(м(1)) — v(м(2)) = V(м(1) — м(2)) является решением задачи (5) при Ш =0 и и = м(1) — м(2). Тогда по неравенству (27)

(и(1) — «(2))|к(п) < и-1М\и(1) — и(2) |Ь2(П). (28)

Аналогичным образом, учитывая оценки (18), (28), получаем неравенство

^'(и«) — Vw(и(2))Нi2(п) < divv(M(1)) — divv(u(2))||L2(п) <

< п (и(1) — И(2))||х2(п) < п3/2Ци(1) — и^п).

(29)

Из двух последних оценок и неравенства (17) следует, что в соотношениях (26) можно перейти к пределу при п ^ го. В результате получаем, что

,]2(и + к) — Ми) = (w(u), ^^(п) + 1 Иу ^^(к)|^2(п) =

2

(— div '(и),к)Ь2(п) + 2 ИУ ^^ (к)||L2(П),

причем, с учетом оценок (18), (28)

Нdiv^^(к)|к(п) < п\\ЩЬ2(п).

Условие Липшица для 32 (и) устанавливается с применением неравенства (29):

||^(и(1)) — ^(«(1))|ь2(п) = |^У'(И(1) — ^^(п) <

< Vw(и(1) — и(2) )|х2(п) < ^2|«(1) — «(2)|к(п).

Таким образом, справедлива следующая

Теорема 3. Пусть выполнены условия теоремы 2; тогда функционал 32(и) дифференцируем по Фреше на и2, его градиент определяется по формуле (9), удовлетворяет условию Липшица с константой Ь2 = и-2 п2.

5. Модифицированный градиентный метод наискорейшего спуска. Существуют различные способы выбора величины а^ (см. [6], [7]). Наиболее быструю сходимость дает метод наискорейшего спуска, в котором а^ определяется из условия

(ак)= тгп ¡к(а), Д (а) = Зг [Рщ(ик — а/г (ик))) , I = 0,1, 2. а>о \ /

(30)

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

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

Если ик € и1, то Рц1 (ик) = ик, поэтому, с учетом обозначения с1к = Ри1 (3 (щ)), имеем Д (а) = 3 (ик -a.dk) = Н^у v(ufc) - а ^^ (4 )|Ц2(П)/2 =

= v(uк)||2(п)/2 - аv(ик), v((1к))ь2(п) + а2|а1уv((1к)||2(п)/2, откуда точкой минимума является

, (а1уv(uk), л^( ¿к))й2(п)

ак =-5- , (31)

1 1 ^(4) 11 ^(П)

а (6) принимает вид

ик+1 = ик - акАк. (32)

Применение метода наискорейшего спуска для рассматриваемых задач встречает затруднение, связанное с тем, что функционалы ^(и) не удовлетворяют условию ограниченности множества Лебега М^с = {и € ^ : 3(и) < с}, которое используется при доказательстве теорем сходимости метода наискорейшего спуска. Оказалось, что эту трудность можно преодолеть, если а к выбирать по формуле

ак = тгп

ак ,1

(33)

где 7 - параметр метода, а ак определяется как в методе наискорейшего спуска по формуле (31).

Поскольку предлагаемый метод может быть использован и в других задачах оптимизации, сформулируем утверждение в виде теоремы в абстрактном гильбертовом пространстве Н. Введем обозначения: J* = inf J (и), U С Н, U* = {и Е U : J (и) = J*}, С1,1 (U) -

и

множество дифференцируемых функционалов, градиент которых удовлетворяет условию Липшица.

Теорема 4. 1 Пусть U - выпуклое, замкнутое множество из гильбертового пространства Н с нормой |1 ■ ||, J (и) Е С1,1 (U) - выпуклый функционал (J* > —ю), множество U* непусто и ограничено, последовательность {икопределена по формулам (33), (6)2 и выполнены условия:

<х

El I J' ) 11 2 ^ Ь1, (34)

к=0

0 <ак < 62; (35)

тогда последовательность {ик}^0 минимизирует функцию J (и) на U и слабо в Н сходится к множеству U*.

Доказательство. Обозначим p(u,,U*) = min ||и — v||; тогда по определению оператора

veu*

проектирования

р2 (uk+1, U*) = ||ий+1 — Ри* (и^)||2 ^ Цик+1 — Ри* (ик)||2

1 Данная теорема с доказательством приведена в работе [10]. Однако утверждение о том, что в случае, если область и является подпространством или всем пространством, то условие (34) можно отбросить, доказано для конкретного функционала. В данной работе это установлено для любого выпуклого функционала 7(и) € С 1'1(и).

2При отсутствии в обозначении функционала J(и) нижнего индекса, аналогичное предполагается и в формуле (6).

= \\Pu (Uk - акJ' (ик)) - Ри (Ри* (uk))\\2 ^ \\uk - акJ' (ик) - Ри* (ик)\\2 =

= р2 (uk, U*) + ак \\J' (uk)\\2 - 2ак (J' (uk) ,Uk - Pu* (uk)). (36)

Воспользовавшись необходимым и достаточным условием выпуклости дифференцируемого функционала на выпуклом множестве U [6] J (u) - J (v) > (J' (v) ,и - v) У и, v G U, полагая в нем v = Uk, и = Pjt (щ), получаем

(J' (ик), Uk - Pu* (uk)) > J (uk) - J (Pu* (uk)) = J (uk) - J* > 0, (37)

откуда неравенство (36) принимает вид

p2 (uk+i, U*) - p2 (uk, U*) ^ a2k \\J' (uk )\\2 .

Суммируя последнее неравенство от 0 до т - 1 (т > 0), и учитывая условие (34), получаем

т— 1

Р2 (Um, и*) ^ ^ а2к \\J' (Uk )\\2 + р2 (щ,и*) ^ 62 bi + Р2 (щ,и*). к=0

Таким образом, последовательность |мк}^=0 ограничена в Н, а из условия (34) следует, что lim, \\J' (ик)\\ = 0; тогда из неравенства (37) следует, что последовательность (ик}^=0

к^-те

минимизирует функционал J (и). Таким образом, последовательность (и к}те=0 - ограниченная и минимизирующая J (и) на U.

Обозначим через W множество выпуклых комбинаций последовательности (ик}^=0, то есть множество точек и, представимых в виде:

те те

и = Y1 акUk, ак > 0, ак = 1. к=0 к=0

Используя [7, гл. 4, §1, теорема 5], легко показать, что W С U и поскольку U - замкнутое множество, замыкание W множества W также принадлежит U.

Последовательность (ик}те=0 минимизирует функцию J (и) на U и, следовательно, минимизирует J (и) на W. Из доказанного следует, что J* {W^j = inf J (и) = J* = inf J (и),

ueW ueu

W* = {u G W : J (u) = J*} G U*. Из ограниченности последовательности (uk}^=0 следует ограниченность множества W. Согласно [6, гл. 1, §3, теорема 6] выпуклый, полуограниченный снизу функционал J (и) на ограниченном, выпуклом, замкнутом множестве U из рефлексивного банахового пространства имеет непустое множество точек минимума U*, и любая минимизирующая последовательность (ик}те=0 слабо сходится к U*. Из слабой сходимости последовательности (ик}^=0 к W* следует ее слабая сходимость к U*, теорема доказана. □

Замечание 1. Если множество W компактно, то имеет место сильная сходимость. Здесь можно воспользоваться [6, гл. 1 §3, теорема 1].

Замечание 2. Если U - подпространство гильбертового пространства Н, Pjj - оператор ортогонального проектирования на это подпространство, то Uk+i = Uk - Pu J' (uk). В этом случае соотношение (36) можно записать в виде:

р2 (uk+i, U*) = р2 (ик, U*) + с?к \\PuJ' (ик)\\2 - 2а,к (PuJ' (щ) ,Uk - Pu* Ю).

Учитывая, что (PuJ' (ик) ,ик - Ри* (ик)) = (J' (ик) ,ик - Ри* (ик)), легко видеть, что утверждения теоремы справедливы, если вместо условия (34) выполняется условие

те

^ \\PuJ' (Uk)\\2 <bi. (38)

к=0

Теорема 5. Пусть и - подпространство или все гильбертово пространство Н, 3(и) € С 1,1 (и) - выпуклый функционал, множество и* непусто и ограничено, последовательность [ик}£=0 определена по формулам (33), (6), тогда последовательность [ик} минимизирует функционал 3(и) на и и слабо в Н сходится к множеству и*.

Доказательство. Покажем, что если 3(и) € С1,1 (и), где и либо подпространство, либо все пространство Н и параметр ак определяется по формулам (33), (30), то выполнено условие (38) и, следовательно, утверждения теоремы верны без условия (34). Для этого воспользуемся известным неравенством, справедливым для функций из С1,1 (и) (см. [7, гл. 2, §3, лемма 1]):

| 3 (и) - 3 (г;) - (3' (г;) ,и - г;)| ^ Ь \\и -V||2 /2 Уи,у € и,

где Ь - константа Липшица для градиента 3 (и) функционала 3 (и). Полагая в нем V = ик, и = и^+1 = ик - аРи3' (ик), получим

3 (ик) - 3 (<+1) = 3 (ик) - 3 (ик - а,Ри3' (ик)) >

> а (/ (ик), Ри3' (ик)) - а2Ь |Р^3' (ик)||2 /2.

(39)

Учитывая, что оператор Ри - оператор ортогонального проектирования на подпространство, получаем, что (3 (ик), Ри3 (ик)) = \\Ри3' (ик)\\, а из неравенства (39) следует:

3 (ик) - 3 (<+1) > а (1 - аЬ/2) \\Ри3' (ик)\\2 . (40)

Полагая а = 1/Ь, получаем

3 (ик) -3 «+0 > \\Ри3' (ик)\\2 /(2Ь).

Предположим, что ак ^7, тогда ак = ак и поэтому при а = 1/Ь справедливы неравенства

2

3 (ик) - 3 (ик+1) > 3 (ик) - 3 «+1) > Ри3 (ик) /(2Ь). (41)

Предположим теперь, что ак > 7, тогда ак = 7. Рассмотрим два случая: 7 > 1/Ь и 7 < 1/Ь. Учитывая, что на интервале (0,а^ функция Д (а) убывает, в первом случае вновь получаем неравенства (41). Во втором случае (7 < 1/ Ь)

7(1 -7Ь/2) > 7/2.

Таким образом, учитывая неравенства (40), (41), в любом случае получаем оценку

2

3 (ик) - 3 (ик+1) > с Ри3 (ик) ,

где с = тгп [7/2, 1/(2Ь)].

Из последней оценки следует, что последовательность [ 3 (ик )}^0 монотонно убывает, ^ / 2

ряд £ | | Ри3 (ик)|| сходится и имеет место искомая оценка

к=0

оо

2

| Ри3 (ик) ^ с-1 (3 (ик) -3*).

3=к

Таким образом, выполнено условие (38), и теорема доказана. □

6. Итерационные процессы. Теорема сходимости. С учетом теоремы 2 далее будем говорить только о задачах I = 1, 2. Покажем, что функционалы 3\(и) удовлетворяют всем условиям теоремы 5.

Нетрудно убедиться, что функционалы 3 I выпуклы. Действительно при У а Е [0,1] и

и,ь Е

3[(аи + (1 — а)ь) = ||а + (1 — а) у(^) |||2(П)

= а2||&уУ(И)||2(П) + (1 — а?ИГу^^(П) + М1 — а) (аГуv(u), аГу^^(П) = а^гуу(и)^^(П) + (1 — аОИ^^^(П) — а(1 — аОИ^у(и) — ^^(П) ^

^ аЗ\(и) + (1 — а)31 (у).

Выполнение условия принадлежности функционалов 3\ классу С1,1 (и) следует из результатов разделов 2, 3, 3%* = т/ 3%(и) = 0 > — го, множества и^* = {и Е III : З^и) = 3%,*}

и

состоят из одной точки и* = р, где р есть искомое давление.

Из теоремы 5, таким образом, следует, что последовательности {ик}£=0, определенные соотношениями (6), (31), (33) при I = 1, 2 слабо сходятся в соответствующих пространствах Н^П) и Ь2(П) при любом начальном приближении.

Теорема 6. Пусть f Е Х2(П), 5 Е С2; тогда последовательность {ик}^=0, определенная формулами (6), (33), где I = 1, 2, при любом начальном приближении и0 Е сходится к и* = р слабо в Н1 (П) и при I = 1 - сильно в Ь2(П). Последовательности {ь(ик)}^0 при I = 1, 2 сильно в Н(П) сходятся к V, где р и V - решение задачи (1)-(3).

Доказательство. Как было показано, первое утверждение сразу следует из теоремы 5. Справедливость второго утверждения вытекает из следующих соображений. Известно (см., например, [1, гл. 1, §1, лемма 8]), что если область П ограничена и ее граница кусочно-гладкая, то Н 1(П) вкладывается компактно в Ь2(П) (т.е. ограниченное в Н 1(П) множество компактно в Ь2(П)). Из компактности и слабой сходимости в Н (П) последовательности {^(ик)}д= следует ее сильная сходимость в Х2(П).

Заметим далее, что разность — V, где = ч(ик), а V - решение задачи (1)-(3), является решением задачи

—и А (V к — V) = —V (ик — р) , (V к — V) ^ = 0.

Умножая последнее уравнение на V к — V и интегрируя по частям, получаем

№(Ук — у)||£2(п) = — р-1^(ик — р), V к — у)£2(п) = "-1(ик — Р, Ук)ь2(П). (42)

Последовательность {ик}<к=0 ограничена в Ь2(П) и минимизирующая функционал 3\(и), поэтому ^у V к ^ 0 при к ^ го. Откуда следует, что (ик — р, Ук)ь2(П) ^ 0 при к ^ го,

поэтому из соотношений (42) следует справедливость последнего утверждения теоремы.

□

Установленная теорема сходимости позволяет сформулировать итерационные процессы в дифференциальной форме для решения задачи (1)-(3). Алгоритм для 31(и)

1. Выбрать начальное и0 Е и1, например, и0 = 0.

2. Найти скорость Ук (ик) как решение векторной задачи Дирихле (5).

3. Найти сопряженное состояние Wк(ик) из решения векторной задачи Дирихле (10).

4. Определить Зк = Ри1 (р('к)), где р('к) - решение скалярной задачи Неймана (25).

5. Определить Vк(Зк) из решения векторной задачи Дирихле (5) при { = 0.

6. Вычислить ак по формулам (31), (33) для найденных Ук(ик), Ук(Зк).

7. Пересчитать управление ик+1 по формуле (32) для найденных ак, Зк и перейти к шагу 2.

Алгоритм для 32(и) идентичен алгоритму 31(и) за исключением шага 4:

4 = Рщ (- Шк(ик)).

Операторы проектирования для обоих алгоритмов вычисляются по формуле

Рщ (и)=и - (щ 1)мп) Уи € Щ, / = 1, 2.

(43)

7. Комбинированный градиентный метод. Состоит в том, что попеременно делается несколько шагов первого градиентного метода (т.е. 3'(и) = 31 (и)), а затем несколько шагов второго (т.е. 3'(и) = 32(и)). Расчеты ниже показывают, что комбинирование методов обеспечивает ускорение сходимости. Этот факт имеет простое объяснение. Функционалы 31 (и), 32(и) достигают минимума через обращение в нуль в точке и* = р тогда и только тогда, когда сопряженное состояние ш(и*) = 0. Действительно, если ш(ик) = 0, то в разложении ш(ик) на градиентную и соленоидальную составляющие, по крайней мере, хотя бы одна часть не равна нулю. Пусть, например, не равна нулю градиентная часть ш(ик), тогда спуск по градиенту 31(ик) обеспечит уменьшение функционала. Комбинированный метод состоит в последовательном переходе от использования одного градиента на использование другого, как только величина \\3/(ик)\\ ( 1 = 1, 2) станет меньше заданной величины.

Заметим, что области определений функционалов и (1 = 1, 2) различные, поэтому требуется убедиться, что по окончании использования одного метода градиентного спуска мы получаем управление ик, принадлежащее области определения другого функционала. Это условие будет выполнено, если начальное приближение и0 € и1. Чтобы в этом убедиться, достаточно показать, что если и0 € Н 1(П), то 3'2(и0) € Н 1(П). Действительно, если и0 € Н 1(П), то при условии гладкости границы Б € С2 в уравнении (5) правая часть f - Уи0 € Х2(П), поэтому его решение у(и0) € Н (П). Откуда следует, что правая часть уравнения (10) принадлежит Х2(П), поэтому его решение ш(и0) € Н2(П) и, следовательно, 32(и0) = - ^у ш(и0) € Н (П). Таким образом, все члены последовательности ик, полученные комбинированным методом, содержатся в и2.

Если применение градиента 32 (и) заканчивается, начиная с некоторого номера итерационного процесса, то, в соответствии с теоремой 6, итерационный процесс сходится слабо в Н 1(П) и сильно в Ь2 (П). В противном случае можно доказать слабую сходимость итерационного процесса в £2(П). Последовательности векторов скорости [\(ик)}£=0 в силу теоремы 6 сильно сходятся в пространстве Н^(П).

8. О задаче с неоднородными граничными условиями. Приведем некоторые замечания о возможном ослаблении условий на границу и применении разработанных методов к решению задачи Стокса с неоднородными граничными условиями.

Условие на гладкость границы области Б € С2 используется лишь для применения второго энергетического неравенства (12). Однако, как отмечается в [8], условие на границу можно значительно ослабить, и этим ослабленным условиям удовлетворяют, например, любой многогранник или произвольная выпуклая область.

Для решения задачи Стокса с неоднородными граничными условиями

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

1/Ду = Г-Ур, = 0,

V = 0,

(Р, 1)Ь2(П) = °

(44)

(45)

(46)

управлению и. Если управлению и дается приращение к, то приращение вектора скорости v(Vh) = v(V(u + к)) - v(Vи) определяется решением однородной краевой задачи как и в ранее рассмотренных случаях. Теперь нетрудно видеть, что выкладки, приведенные при доказательстве теорем 1, 2, 3, остаются без изменений.

Таким образом, функционалы 31 (I = 1, 2), рассматриваемые в пространствах Н 1(П) и Ь2(П) и в случае неоднородных краевых условий, дифференцируемы и удовлетворяют условию Липшица, кроме того, сохраняется формула (31) для определения параметров ®к.

Для обоснования сходимости последовательности {ик}^=0, построенных указанными выше методами, необходимо показать, что множество точек минимума функционалов 3\ непусто и ограничено. Для этого достаточно убедиться, что решение задачи (44)-(46) существует.

Существование решения задачи (44)-(46) в классе V € Н2(П), Vp € ¿2(П) при условии, что область П ограничена, f € Х2(П), Б € С2, 0^ € Н1+2, (<0,п)ь2(вГ) = 0 гарантирует [1, гл. III, §5, теорема 3]. Следовательно, при указанных условиях, справедливы утверждения теоремы 6.

9. Конечно-разностная аппроксимация. Для построения сеточных аналогов дифференциальных итерационных процессов из раздела 5 недостаточно воспользоваться произвольными аппроксимациями операторов ^у и V и гарантировать при этом сходимость, с чем и столкнулись авторы при численном моделировании. Однако, если сначала аппроксимировать исходную задачу (1)-(3) в конечномерных пространствах Соболева и построить для этой задачи градиентный метод аналогично рассмотренному выше дифференциальному случаю, удается получить сходящийся сеточный итерационный процесс, для которого справедлива теорема 6.

Как и в работе [11] будем рассматривать двумерную прямоугольную область Пн = (0, /1) х (0,12) и узловую равномерную неразнесенную сетку = {xí,j = (ък1,]к2), г=1,..., Ы1, ] = 1, . . . , К2} с границей Бн (состоит из множества граничных узлов, за исключением угловых точек). Здесь ка = Iа/(Ыа + 1) - равномерный шаг, а К - число внутренних узлов сетки по направлению а = 1, 2. В дальнейшем ограничимся случаем равномерной сетки N = К2 = N в квадрате /1 = 12 = I (к1 = к2 = к).

Определим скалярное произведение в конечномерном пространстве Ь2(Пн):

к2 N N

(у, ^ Ь2Ш = к2 Е у™ + Т Е У™ = к2 Е Е Уч

к2

2

Хг^ еЯи í=1 3=1

(N N N N \

Уí,N+1 +1 + Е У0,з20,] + Е УN+1,jzN+1, з I

í=1 í=1 3 = 1 3 = 1 )

У%,0 zí,0 У^+1 +1

где у^ = у(х^э), г^у = г(х^3).

Под согласованной аппроксимацией операторов ^Ун и Vh понимается сохранение дифференциального свойства интегрирования по частям:

(VhP, ^(Пн) = -(Р, (Пн) Vv € Н01(Пн). (47)

Соотношению (47) удовлетворяют направленные разности первого порядка [11]:

(1) (1) (2) (2)

diУhVгз = ^ , í-lJ + íJ , íJ-1, г = 1 ...К + 1, = 1 ...К + 1, к к

VhPг,з = (^+1"к- 1КЗ , ^ ) , . = 0 ...К, 3=0 ...К.

Очевидно, что в этом случае оператор -Дн = - divhVh есть стандартный пятиточечный разностный шаблон для задачи Дирихле [12].

Таким образом, сеточным аналогом задачи (1)-(3) является задача

— иАну = ^ — Vhp, у|

Я

0,

(48)

&уку = 0, (49)

(Р, 1)ь2(п,) = 0, (50)

для которой выпишем итерационные процессы из раздела 5 в удобном при программной реализации виде. При этом вместо 31(и) будем теперь рассматривать 30(и). Алгоритм для /0(м)

— уАьук = 4 — ик, V к

вн

0;

—иА^к = к, w к к = 0;

дрк

— Арк = — ^у^ к,

дп

(wk ■ п^вн;

(51)

вн

к = Ри1 ( к ); —рАъук = —V},Зк, Vк|

К

0;

ак = тгп

к Л

а1

(¿Гу^Ук, ¿Гу^У к)ь2(пн),

||СГу}Ук ^2(Пн)

«к+1 = «к — «к Vhdk.

Алгоритм для 32(и)

0;

—уАьук = 4 — Vuk, V к |вн — иА^к = ^¿ГУ^ук, wк |вн = 0; Рк = — ¿Гу^ к; Зк = ри1 (Рк);

—рАъУк = —V} Зк, Ук |вн = 0;

ак = тгп

«к Л

ак =

(¿Гу^Ук, сИу^к)ь2(пн),

||сГу}ук ^2(Пн)

«к+1 = ик — а.кЗк.

Для обоих алгоритмов

риг (Рк) = Рк —

(рк, 1)Ь2(Пн) (1, 1)Ь2(Пн )

где, как нетрудно убедиться, (1,1)ь2(пн) = кк.

10. Решение вырожденной задачи Неймана. Хорошо известно, что условием разрешимости дифференциальной задачи Неймана вида

—Ар = ! д£.

' дп

является равенство

(Л 1)Ь2(П) + (д, 1)Ь2(5) = 0, которое для задачи (24) принимает вид формулы Остроградского-Гаусса

(¿гу wk, 1)Ь2(П) = (^к ■ n), 1)Ь2(в).

Справедливость последней формулы для сеточного случая проверяется непосредственно, если скалярное произведение в Ь2(Зи) определено как

N N

{д, ^(Я) = Ь ^ (90,% + д?,г ) + Ь ^ 1'3 + 93, А +

г=1 з=\

Ь

+ ^^о,о + 5*1,0 + 52,N+1 + да, N+1 + 52,0 + д1, N+1 + до, N+1 + #з,о) Уд е ^(^и).

Введем обозначение компонент вектора w к = (и^, и^). По определению дрк/дЩ3 = (VнРк ' ^)1я1г, поэтому граничное условие для задачи (51) определяется из равенства (VиРк — Wk)|я = 0:

(Рк)ц — (Рк)о _ Ь = (Рк ) г,1 — (Рк ) г,0

3 = (,„(1)) (Рк)N+2^ - (Рк)N+Ц = (7 (1)) .

" = (тк к^ -Т- = (тк )N+l,j, 3 е М*.

Ь

(Рк) ¿,N+2 — (Рк )

/ (2)ч (Рк)¿,N+2 - (Рк)¿,N+1 , (2) \ ■

К ))^ -т-= К ))г^+ъ ге 1,М.

Ь ^ Ь у к

Из приведенных формул видно, что требуется введение фиктивных слоев точек по каждому из направлений с индексами N + 2 за пределами сеточной области ши и , после чего можно записать систему ( N + 2) х ( N + 2) - линейных алгебраических уравнений для решения задачи (51):

Аг = Ь, (53)

где гСЕьь(г,з) = (Рк)а оператор СЕЬЬ(ъ, з) = г^ + 2) + з (г,3 е 0, N + 1) задает построчный формат хранения двумерного вектора в одномерном массиве.

Все элементы строки г от = СЕЬЬ(г, з) (г,3 е 0^ + 1) матрицы А равны нулю, за исключением:

(Ч ЬЗ е 1, N;

А

гот, СЕЬЬ(г, з)

3, г = 0 и ^ е 1, N или г = N +1 и] е 1, N;

3, j = 0 иг е 1, N или з = N + 1 иг е 1, N; 2, г = 0 и э = 0 или г = 0 иj = N + 1; к 2, г = N + 1 и з = 0 или г = N +1 из = N + 1.

Агот,СЕЬЬ(г-1, з) = -1, 2 > 1.

А

■гот,СЕЬЬ(г+1, з)

,) = -1, г< N;

Агот,СЕЫ(г, з-1) = -1, 3 > 1; Агот,СЕЬЬ(г, з+1) = —1, 3 ^

Соответственно, вектор правой части вычисляется по формулам

0, г е 1, N, 3 е 1, N; (w<k1))о,з, г

1, 3 е 0, N + 1;

Кот = ) 1,з + Ь{ (тк2))¿,N+1, = N +1, г е 0, N +1;

'к1):

(2)

+и, г = N +1, з е 0, N +1;

. МЛг,о, 3 = 1, ге 0, N +1.

Система (53) имеет разреженную структуру и может быть эффективна разрешена подходящим итерационным методом с предобусловливанием.

11. Численные эксперименты. Традиционно, первичная верификация численных методов осуществляется на задачах с известным аналитическим решением. Применительно к задаче (1) - (3) это означает, что задавая соленоидальное V* и р* из условия (р*, 1)ь2(п)

мы можем вычислить сеточную правую часть fh = (—uAv* + Vp*)h и найти вектора , Uk в результате сходимости алгоритмов из раздела 1.

Рассмотрим тестовый пример из работы [4]. В области П = \-ж/2, 3^/2] х [—ж,ж\ решается задача Стокса (1)—(3), где для v = 1 задается соленоидальное аналитическое решение v* (х, у) = {(1 + sinx)sin у, cosx(1 + cos у)}, р*(х, у) = sinx cos(2y), для которого v*|s = 0.

Задача (48)-(50) решалась на последовательности сеток 31 х 31, 63 х 63, 127 х 127, 255 х 255 с помощью алгоритмов J2(u), Jo (и), а также их комбинированием. Задачи Дирихле для уравнения Пуассона обращались методом быстрого преобразования Фурье [12], с чем и связан выбор числа внутренних узлов сеток по закону N = N2 = N = 2т — 1 (т = 5, 6, 7, 8). Задача Неймана (53) решалась итерационно методом сопряженных градиентов с диагональным предобусловливанием. Во всех алгоритмах итерации продолжались до тех пор, пока ||pk — Pk-i||c(nh) > 10-6 или ||div^v k||c(oh) > 10-6, где || • ||c(nh) - равномерная сеточная норма, т.е. максимум среди всех модулей функции на сетке ш^. В качестве вектора начальных приближений использовано и0 = 0.

Результаты вычислений с помощью разработанной авторами программы на языке C+—+ приведены в табл. 1, 2; на рис. 1, 2 представлена динамика убывания нормы дивергенции вектора скорости.

Сетка Итерации ||v k — v*Wc (Oh) W rPk— p*!c(Oh) Невязка |f — Vhpk + Ahvk||c(oh)

31 х 31 200 1.92952e-11 3.05475e-5 1.54609e-9

63 х 63 242 7.96903e-12 3.81565e-6 2.72057e-9

127 х 127 262 3.3259e-12 4.76882e-7 5.99933e-9

255 х 255 265 1.41589e-12 5.96158e-8 1.26728e-8

ТАБЛИЦА 1. Результаты расчета по алгоритму J2(u)

Сетка Итерации ||v k — v*Wc (Oh) WPk — P* | c(Oh) Невязка ||fh — VhPk + Ahvk||c(Oh)

31 х 31 167 1.80279e-11 1.7414e-6 1.15983e-9

63 х 63 143 7.09932e-12 1.56985e-5 2.04383e-9

127 х 127 163 4.07434e-12 1.0086e-5 5.90967e-9

255 х 255 219 2.13486e-12 5.4197e-6 1.1914e-8

ТАБЛИЦА 2. Результаты расчета по комбинированному алгоритму, где после одного шага и) итерации продолжаются по алгоритму 32{и)

Расчеты по алгоритму 30(и) не приводятся ввиду обнаруженной вычислительной неустойчивости соответствующего расчета.

Влияние параметра 7 на сходимость модифицированного метода градиентного спуска с вычислением параметра ак по формулам (33), (31) проиллюстрировано на рис. 2. За основу был выбран комбинированный метод, где после одного шага и) дальнейшие итерации производились по алгоритму 32{и). Оптимальным значением 7 на рассмотренных сетках является число в пределах 9-11, которое, скорее всего, будет лежать в других диапазонах при других условиях задачи. Сравнение соответствующей табл. 3 с табл. 2 показывает сокращение числа итераций.

12. Заключение. В заключение отметим, что, во-первых, два построенных итерационных процесса из раздела 1 записаны в дифференциальной форме и инвариантны к размерности рассматриваемой задачи Стокса (п = 2, 3). Во-вторых, редуцирование исходной задачи к серии существенно более простых задач Дирихле и Неймана, позволяет использовать для их численного решения известные эффективные сеточные методы.

0 50 100 150 200 250

-0-31x31 -6-63x63 —0—127x127 -0-255x255

Рис. 1. График изменения log10 (&Уи Vk) для алгоритма 32(и) в зависимости от номера итерации

0 50 100 150 200

Рис. 2. График изменения log10 (^Уи V к) для комбинированного алгоритма в зависимости от номера итерации

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

Алгоритм 32(и) проще при программной реализации, однако сходится медленнее, чем комбинированный метод.

Сравнение с численными результатами из работы [4] показывает, что даже при применении базового метода наискорейшего спуска предложенные в настоящей работе алгоритмы

gamma

31x31 -а-63x63 -°-127x127 -°-255x255

Рис. 3. Зависимость числа итераций комбинированного метода от 7.

Сетка Итерации IIVk - V*||c(Qfe) ||uk- p*h (nh) Невязка ||f - VhUk + AhVk ||c(nfe)

31 х 31 113 2.53668e-12 2.00158e-6 1.56719e-9

63 х 63 121 0.999382e-11 1.00493e-5 2.94942e-9

127 х 127 133 0.497812e-11 0.502867e-5 4.47148e-9

255 х 255 189 0.248806e-11 0.251476e-6 9.65243e-8

ТАБЛИЦА 3. Результаты расчета по комбинированному алгоритму при 7 = 10

затрачивают примерно в два раза меньше итераций для достижения точности по скорости

во втором знаке на сопоставимых сетках.

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

1. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М: Наука, Главная редакция физико-математической литературы, 1970.

2. Темам Р. Уравнения Навье-Стокса. Теория и численный анализ. М.: Мир, 1981.

3. Экланд И., Темам Р. Выпуклый анализ и вариационные проблемы. М.: Мир, 1979.

4. V.l. Agoshkov, C. Bardos, S.N. Bideev Solution of the Stokes problem as an inverse problem: Preprint N9935. Centre Math. Leuers Applic. Cachan: E.N.S. de Cachan, 1999.

5. Пальцев Б.В., Соловьев М.Б., Чечель И.И. О развитии итерационных методов с расщеплением граничных условий решения краевых и начально-краевых задач для линеаризованных и нелинейной систем Навье-Стокса // ЖВМиМФ. 2011. 51:1. C. 74-95.

6. Васильев Ф.П. Методы решения экстремальных задач. М.: Наука, 1981.

7. Васильев Ф.П. Численные методы решения экстремальных задач. М.: Наука, 1988.

8. Ладыженская О.А. Краевые задачи математической физики. М: Наука, Главная редакция физико-математической литературы, 1973.

9. Ладыженская О.А. О связи задачи Стокса и разложений пространств W2 и W2 ^ // Алгебра и анализ. 2001. 13:4. C. 119-133.

10. Голичев И.И. Модифицированный градиентный метод наискорейшего спуска решения линеаризованной задачи для нестационарных уравнений Навье-Стокса // Уфимский математический журнал. 2013. T. 5. № 4. C. 60-76.

11. A.G. Churbanov, A.N. Pavlov, P.N. Vabishchevich Operator-splitting methods for the incompressible Navier-Stokes equations on non-staggered grids. Part 1: first-order schemes // International Jornal for Numerical methods in fluids. 1995. V.21. P. 617-640.

12. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989.

Голичев Иосиф Иосифович Институт математики c ВЦ УНЦ РАН ул. Чернышевского, 112 450008, г. Уфа, Россия

Уфимский филиал Финансового университета при Правительстве РФ ул. Революционная, 169 450005, г. Уфа, Россия E-mail: [email protected]

Шарипов Тимур Рафаилевич, Научно-производственное предприятие «АТП» ул. Кузнецовский Затон, дом 22, корпус 2 450103, г. Уфа, Россия E-mail: [email protected]

Лучникова Наталья Иосифовна,

Уфимский филиал Финансового университета

при Правительстве РФ

ул. Революционная, 169

450005, г. Уфа, Россия

E-mail: [email protected]