Метод решения задачи Дирихле для эллиптического уравнения 2-го порядка в полигональной области на последовательности адаптивно измельчаемых сеток Текст научной статьи по специальности «Математика»

УДК 519.632.4

Е. С. Николаев, О. В. Шишкина

МЕТОД РЕШЕНИЯ ЗАДАЧИ ДИРИХЛЕ ДЛЯ ЭЛЛИПТИЧЕСКОГО УРАВНЕНИЯ 2-го ПОРЯДКА В ПОЛИГОНАЛЬНОЙ ОБЛАСТИ НА ПОСЛЕДОВАТЕЛЬНОСТИ АДАПТИВНО ИЗМЕЛЬЧАЕМЫХ СЕТОК

(лабораторияразностных методов факультета ВМиК)

1. Введение. Методы конечных элементов (МКЭ) решения краевых задач для дифференциальных уравнений, использующие адаптивные стратегии измельчения сетки и локального изменения порядка элементов (h- и р-версии), активно развиваются с конца 70-х годов. Эти методы основаны либо на прямом использовании локальных апостериорных оценок погрешности, математический аппарат получения которых был развит в работах [1-6], либо на вычислении индикатора коррекции — величины, оценивающей степень изменения выбранной характеристики приближенного решения при добавлении иерархическим образом одной или нескольких новых степеней свободы [7-11]. Современное состояние исследований по технике построения апостериорных оценок отражено в обзорных статьях [12-15], а различия между этими классами адаптивных методов подробно обсуждаются в работе [8].

Кроме указанных отметим методы, генерирующие адаптивное распределение узлов сетки при заданном их числе (r-версия МКЭ) [16-18], а также основанные на перестройке заданной начальной сетки путем прямой минимизации функционала относительно координат узлов сетки [19-21] методы для задач, допускающих вариационную формулировку.

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

2. Задача минимизации и ее сеточный аналог. Требуется решить задачу минимизации функционала: в полигональной области G = G + dG с границей dG найти функцию и*(х) такую, что

и* = arg min J(u), и* G Н = {и G W-j (G), и(х) = щ(х), х G dG} , (1)

иеН

где J(u), билинейная a(u,v) и линейная b(v) формы определяются следующим образом:

1 С ^ д д С

J(u) = -а(и,и) - b(u), a(u,v)= / ( X] + qUV) dx' b(v) = / fvdx. (2)

1 «,/3=i Xß i

Предполагается, что выполнены условия:

1) матрица k = {kaß(x) ^ симметрична и равномерно по х G G положительно определена, q = q[x) — неотрицательная функция;

2) функция щ непрерывна на dG и допускает продолжение в W^iG), т.е. существует функция v G W\ (G) такая, что ее след на dG совпадает с щ.

В силу сделанных относительно к и q предположений форма а(-, •) симметрична и на функциях, обращающихся в нуль на dG, коэрцитивна.

Если ка/з, q G L00{G)^ а / G L/2(G), то решение задачи (1), (2) существует, единственно и является обобщенным решением задачи Дирихле для эллиптического уравнения 2-го порядка

2 д f Ou \

- ( ка117й~ ) + qu = u(x) = ub(x), х G dG. (3)

a, /3 = 1 Xa ^

Если же kafj G G1 {G) и g, / G C(G), a и* G C2(Gr) П C(Gr), то и* — классическое решение задачи (3).

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

_ "< .

Пусть в области построена триангуляция G = [j Тг, содержащая nt треугольников Тг и nv вер-

¿=1

шин, принадлежащих множеству узлов U = из + доз, где

UJ

= {Xi G G, 1 ^ i ^ га}, диз = {Xi G dG, 1 ^ i ^ m}, nv = n + m.

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

Пусть {v:'j}J_1, {фj}j~i — наборы кусочно-линейных непрерывных базисных функций, ассоциируемых с узлами из из и доз. Функции <~Pj{x) и ipj(x) равны единице в узлах, с которыми они ассоцииру-

° Г " 1

ются, и нулю в остальных узлах. Обозначим через Нп = < w = Wj'fj \ множество кусочно-линейных _ i=1

непрерывных на G функций, обращающихся в нуль на dG.

Соответствующая (1) конечноэлементная задача минимизации формулируется следующим образом: в области G найти функцию и(х) такую, что

о га

и = arg min J(u), и G НПгГП = <u = w + g, -u? G Hn, 9 =/ иь{х])Ф] \ ■

3 = 1

Известно, что функция u(x) находится из условия Галеркина

о

a(u,v) = b(v) Vv G Я„. (4)

о ^

Так как базис в конечномерном линейном пространстве Нп образуют функции {Vi}^) то это условие может быть записано в виде системы линейных алгебраических уравнений

Aw = b, (5)

A={a'j}"J=i' w = (wi,...,wn)T, b = (&i,..., bn)T,

a,ij = a(ipj,ipi), bi = b((pi) - a(g, (pi), 1

о

В силу определения <~pi{x), симметрии и коэрцитивности формы а(-, •) на функциях из НП1 а также справедливого для любых векторов w = ..., wn)T, v = (t>i,..., vn)T равенства

n n

Eo __о

WjLpj G Hni V = 2_,vjLPj £ Hn 3=1 3=1

матрица А является разреженной симметричной и положительно определенной. Значение функционала на функции и G НП:ГП1 удовлетворяющей (4), можно вычислить по формуле

п

J(u) = -0,5j2bJwJ+J(g). (6)

3 = 1

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

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

3. Стратегии адаптивного измельчения сетки. Предлагаемые стратегии основаны на использовании двух неотрицательных кусочно-постоянных функций I и Я. Первая из них — индикатор коррекции, значение которой Р на треугольнике Тг оценивает изменение выбранной характеристики приближенного решения, получаемого при размещении на Тг пробного узла и введении в рассмотрение связанной с ним кусочно-линейной функции с носителем, сосредоточенным на этом элементе (как в иерархической /г-версии МКЭ). Вторая функция оценивает точность приближения функции щ(х) кусочно-линейной функцией д(х) на каждом сегменте границы (отрезке между соседними граничными узлами). Значение № на сегменте определим таким образом:

кг =-[Ы~д)2<1х. (7)

тез(З^) .)

ас,-

Стратегия измельчения сетки. Правило отбора треугольников следующее: Тг включается в список Ыв^ отобранных треугольников, если выполняется условие

I* > /Ьаг, Да, = ^ ШЯХ Г + ^

г = 1

где в 6 [—1,1] — параметр, управляющий степенью измельчения сетки. При в = — 1 будут отбираться все треугольники, а при в = 1 — только те, которые имеют максимальное значение индикатора коррекции.

Правило отбора сегментов границы определим следующим образом: дGi включается в список Ыв^ отобранных сегментов, если выполняется условие

Дг ^ тах(дДтах,4), Дтах = тах Д',

где ¡1 6 [0,1] — параметр, с помощью которого можно управлять степенью измельчения сетки, а еи — требуемая точность решения задачи.

В статье рассматриваются два типа индикаторов коррекции — Д и Д. Соответствующие им стратегии измельчения сетки обозначаются через 51 и Индикатор первого типа оценивает уменьшение значения функционала, а второго — изменение решения в добавленном узле. При вычислении индикаторов будем различать внутренние и приграничные треугольники, т.е. такие, одна или две стороны которых (называемые граничными) принадлежат дС. Для каждого треугольника Тг значения индикаторов 1{ и Ц вычисляются при добавлении к множеству функции <£г(х), ассоциируемой с центром этого треугольника. Для приграничных треугольников дополнительно такие же величины вычисляются и при добавлении к множеству •_1 функции фг(х), ассоциируемой с серединой граничной стороны.

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

триангуляцию области. Пополним множество кусочно-линейной базисной функцией срг(х),

ассоциируемой с узлом хг, и определим множества функций

п о т

К+1 = = + К+1,ш = {u = wг+g, wl Е Нгп+1, д = ^Мхз)Фз)-

3=1 3=1

Пусть и 6 Нп,т — приближенное решение задачи (1) на сетке с п внутренними и т граничными узлами, т.е. функция, удовлетворяющая условию (4). Пусть иг 6 Нгп+1 т — приближенное решение задачи (1) на пробной сетке, т.е. функция, удовлетворяющая условию

а(ч\уг)=Ъ(уг) УугЕНгп+1. (8)

Обозначим через и .12 величины

,Гг = J(u) - .УЫ1), ,Г2 = тах Iи(х) - и1(х) I ,

хет<

характеризующие изменение функционала </(и) и решения в треугольнике Т1 при переходе от минимизации функционала на множестве НП)ГП к минимизации его на множестве Нгп+1 т. Теорема 1. Имеют место представления

А = 0,5<+1[%г) - а(и,рг)], Г2 = |<+1| тах|уг(ж) - рг(х)\ , (9)

и- и1 = 1лгп+1 {у1 - р1 , 1лгп+1 = 7 ' -у 7 ' 10

а{(рг, (рг) - а(уг, (рг)

О

где функция уг £ Нп определяется условием

а(уг1у) = а(рг1у) Уг;бЯ„. (11)

Доказательство. Покажем сначала, что для гг = и — и1 имеют место равенства

а{г\р1) = а(и,1рг) - Ъ(рг), ф» = 0 Уг; 6 #„• (12)

о о

Действительно, используя (8) с V1 = (рг, получим первое равенство. Так как Нп С то второе

равенство следует непосредственно из (4) и (8).

о

Докажем справедливость формул (10). Для V = г1 — гигп+1(уг — срг) 6 Нп с учетом (11) и второго из

равенств (12) получим а (и, и) = а(гг,ю) — тгп+1а(уг — = 0. В силу положительности квадратич-

ен

ной формы на отличных от нуля функциях из Нп это равенство имеет место тогда и только тогда, когда V = 0. Отсюда следует первая из формул (10). Подставляя в первое из равенств (12) полученное для гг представление, будем иметь формулу (10) для гигп+1. Покажем, что знаменатель в этой формуле отличен от нуля. Используя условие (11) с V = у\ получим а(срг — уг,рг) = а(срг — уг,рг — уг).

о о

Так как квадратичная форма положительна на отличных от нуля функциях из а <~рг ^ НП1 то

а{р\рг) - а(у,рг) > 0.

о

Используя симметрию формы а(и, и) и соотношение (8) с V1 = г1 6 будем иметь

2,1{ = а(и, и) - а(иг, и1) - 2Ъ(г1) = а(г1, г1) = а(г1, у) - югп+1а(гг, срг),

о

где V = гг + тгп+1срг 6 Нп. Отсюда и из (12) следует формула (9) для Теорема доказана.

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

р = [Ь(р') — а(и, у')]2 р = |%*)-а(И,уг)| 1 2 а((р\(р1) ' 2 а((р\(р1) ' ^ '

Обоснование такого выбора индикаторов коррекции будет дано в пункте 3.3.

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

3.2. Индикаторы коррекции для приграничных треугольников. Обозначим через 1[ и Ц вычисленные по формулам (13) значения индикаторов коррекции для приграничного треугольника Тг. Дополнительно найдем аналогичные величины 1[ и Ц, добавляя новый узел на середину граничной стороны. Эти величины будут характеризовать изменение функционала и решения при замене на этой стороне линейной интерполяции функции иъ при помощи функции д на более точную, кусочно-линейную интерполяцию функцией дг. Значения индикаторов коррекции для приграничного треугольника определим следующим образом:

/* = тах(/',/'), /*= тах(/*,/*).

Если у треугольника Тг две стороны принадлежат границе ЗС, то величины /| и Ц вычисляются для каждой из этих сторон.

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

определим множества функций

п ш

К = = нп,т+1 = {« = ™1 + д\ е Нгп, дг = ^2Мх1)Ф1 + 9гт+1Фг

3=1 3=1

Так как равенство дг(х) = щ(х) должно выполняться для х = хг, то

У'т+1 = Ч(х') - »"*> + *'Ч (14)

Пусть и Е НП)ГП — приближенное решение задачи (1) на сетке с п внутренними и т граничными узлами, т.е. функция, удовлетворяющая условию (4). Пусть иг 6 Нгп т+1 — приближенное решение задачи (1) на пробной сетке, т.е. функция, удовлетворяющая условию

а(и» = Ь(и) Vг;eЯ;. (15)

Как и раньше, обозначим .!{ = -У{и) — .1(и1) ж 3\ = тах \и — иг(х) .

хет< 1 1

Теорема 2. Имеют место представления

А = 91+1 [КФг) ~ Ф1] ~ 0,5(^+1)2 [а(Ф\ Фг) ~ а(у\фг)], (16)

и-иг =дгт+1{уг-фг), 4 = \дгт+1\тах\уг(х) - фг(х)\, (17)

хет•

о

где дгт+1 определено в (14), а функция уг 6 Нп определяется условием

а{у\ю) = а{ф\ю) УуЕНп. (18)

о о

Доказательство. Так как Нп = Нгп, то из (4) и (15) для гг = и — иг будем иметь

а{г\у) = 0 Уу Е Нп. (19)

о

Докажем справедливость формул (17). Для V = гг — дгт+1(уг — фг) Е Нп с учетом (18) и (19) получим а(и, и) = а(гг, и) — дгт+1а(уг — фг, и) = 0. В силу положительности квадратичной формы на отличных от

о

нуля функциях из Нп это равенство имеет место тогда и только тогда, когда V = 0. Отсюда следуют

формулы (17). Учитывая симметрию формы a(u,v), будем иметь

2.J[ = а(и, и) - а(иг, иг) - 2b(zl) = 2 [а(и, zl) - b(z1)] - a(z\ zl) =

= 2[a(u,z') - b(z*)] - a(z\v)+ дгт+1а(гг,фг),

о

где v = z' + дгт+1фг £ Hn. Используя (4) с указанной функцией v, равенство (19) и представление для z%, получим формулу (16) для .J\. Теорема доказана. Осталось определить величины

П = Wm+1 [Ь(Ф1) ~ а(и, " 0,5(^+1)2а(^, | , Ц = \д'т+1 \ ,

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

3.3. Обоснование выбора индикаторов коррекции. Определим билинейную форму

а1 (и, v) = У ( X] + 1UV )

Ti а, /3 = 1

^ ti/ s\ ^ ti/

/3

В силу сделанных относительно матрицы к и функции q предположений справедлива

Лемма 1. Для отличной от нуля на Тг функции v £ W\ (Тг) равенство a'(v, v) = 0 имеет место тогда и только тогда, когда q(x) = 0, v(x) = const при х £ Тг.

Обозначим через (pi , 1 ^ j ^ 3, базисные функции, ассоциированные с вершинами треугольника Тг. Определим матрицу С и вектор z:

a'(Vin ViJ a'(Vin Vis) a'(Vii, v')

С = а*{<Рг 2, «'(У.з» Via) Z = v')

a'(Vis> Vii) О* Via) «'(Via» Via) a'(Vis, V*)

а также обычным образом скалярное произведение векторов: (г, V) = гтлг.

Лемма 2. Матрица С симметрична и неотрицательна, а система уравнений

= г (20)

совместна. Матрица С вырождена тогда и только тогда, когда q(x) = 0 при ж £ Тг, и в этом случае подпространство кегС одномерно и состоит из векторов v = (с, с, с)т.

Доказательство. Симметрия С очевидна, а неотрицательность следует из равенства

a\v,v) = (Cv, v),

где v(x) = У] Vjipi 3 = 1

(ж), v = (ub v2, V3)J

(21)

Пусть матрица С вырождена. Тогда существует ненулевой вектор v такой, что Cv = 0, и поэтому в силу (21) имеет место равенство a'(v, v) = 0. В силу леммы 1 это возможно лишь тогда, когда q(x) = О и v(x) = const при х £ Тг. Отсюда следует, что подпространство кегС состоит из векторов указанного в лемме вида.

Пусть теперь q = 0 на Тг. В силу соотношения <~Pil{x) + <~Pi2{x) + <~Pii{x) = 1, справедливого для х £ Т\ получим, что для вектора v = (1,1,1)т имеют место равенства

Cv= (<¿(^,1), a'(Lpi2,1), а\<р{ 3,1))Т = 0, (z, v) = a\l, = 0.

Отсюда следуют вырожденность С и принадлежность z подпространству imС, ортогональному кегС. Так как z £ im С, то система уравнений (20) совместна. Лемма доказана.

о

Лемма 3. Пусть функция уг £ Нп удовлетворяет условию (11). Равенство

аг(у\уг)= 0 (22)

имеет место тогда и только тогда, когда

аг(¥>>) = 0 Vu £ Я„. Если выполнено условие (23), то уг{х) = 0 для х £ G.

Доказательство. Пусть выполнено (23). Так как носитель функции <рг сосредоточен на Тг, то a(ip',v) = a'(ip',v) = 0. Используя следующее из (11) при v = уг соотношение

а(р\уг)=а(у\уг), (24)

о _

получим а(уг,уг) = 0. Так как форма а(-, •) коэрцитивна на функциях из НП1 то уг(х) = 0 для х £ G, и поэтому равенство (22) справедливо.

Пусть выполнено (22). Тогда либо у* = 0 на Тг и, следовательно, а(рг,уг) = аг(рг,уг) = 0, либо в силу леммы 1 q = 0, уг = const на Тг и поэтому

а(рг, у1) = аг(рг, у1) = J qy\l dx = 0. В силу равенства (24) будем иметь а(уг,уг) = 0. Так как форма а(-, •) коэрцитивна на функциях из

о _

Нп, то уг(х) = 0 для х £ G. Отсюда и из (11) следует (23). Лемма доказана.

о

Лемма 4. Если функция уг £ Нп удовлетворяет условию (11), то верны неравенства

(25)

Здесь 0 ^ 7 < 1 — константа, определяемая следующим образом: 7 = 0, если (23) выполнено, иначе

ИуАУ8)]2

Доказательство. Левое неравенство в (25) следует из (24) и коэрцитивности формы а(-,-). Докажем правое неравенство. Если выполнено (23), то а(рг, уг) = аг(рг, уг) = 0 и правое из неравенств (25) имеет место с у = 0. Пусть теперь условие (23) не выполнено. Тогда в силу леммы 3 будем иметь a'(y'iy') Ф 0- Так как 0 < аг(уг,уг) ^ а(уг,уг) и аг(рг,рг) = а(рг, рг), то справедливо неравенство

[а{р\ f)}2 = 7aV, V^KGA У') ^ f)

с указанным в (26) значением 7. Отсюда с учетом соотношения (24) получим искомое неравенство. Неравенство 7^1 следует из неравенства Коши-Буняковского для интегралов, неравенства Коши

о

для сумм и симметричности матрицы к. Так как ip1 ^ НП1 то 7 < 1. Лемма доказана.

В следующей лемме дана оценка для 7 через величины, не зависящие от функции уг и которые могут быть найдены одновременно с вычислением индикатора коррекции. Лемма 5. Для у из неравенства (25) имеет место оценка

iz,w)

7^7= '. \ < 1, аг(рг, рг)

где w — решение системы (20), если q ф 0 на Тг, и какое-либо ее решение, если q = 0.

Доказательство. Если условие (23) выполнено, то z = 0. Поэтому уг = 0 и совпадает со значением у, указанным в лемме 4 для этого случая. Пусть теперь (23) не выполнено. Из леммы 3 следует, что в рассматриваемом случае аг(уг,уг) > 0. Из (20), (21) и равенства a'(ip',v) = (z,v), где з

V(X) = X) Vj(pi. (ж), v = (vi,v2,v3)T, получим 3 = 1

[а'(р\у')}2 ИУ»]2 (Cw, v)2

< max -—--— = max

аг(УгтУг) ^ г):а; («,«)> 0 аг (V, V) (Су, у)

Так как в силу леммы 2 матрица С симметрична и неотрицательна, то имеет место неравенство (Cw,v)2 ^ (Cw,w)(Cv,v). Отсюда из полученной выше оценки и (26) следует формула для уг. Заметим, что в случае вырожденной матрицы С в качестве уу можно взять какое-либо решение системы (20), так как любые два ее решения отличаются на произвольный вектор, принадлежащий кегС, а вектор г ортогонален этому подпространству.

3

Покажем, что уг < 1. По w = построим функцию w(x) = Wj(pi.(x) ф 0. Так как

з=1

(z,w) = a'(ip',w) = a'(w,w), то для уг имеет место представление, получаемое из (26) при замене уг на w. Завершение доказательства леммы очевидно. Из лемм 4 и 5 вытекает следующая теорема.

Теорема 3. Если выполнено условие (23), то индикаторы коррекции 1\ и Ц, определяемые формулами (13), точны, т. е. Ра = Зга, а = 1, 2. Если (23) не выполнено, то справедливы неравенства

р < р < _-_Т1 Т1 < I < _-_т*

1 _ ihi h Ч \w»+i\ \ 1 _

где постоянная уг определена в лемме 5.

Замечание 1. На треугольнике Тг ненулевые значения имеют только три базисные функции <рч. .

о

Поэтому условие (23) эквивалентно требованию аг {<рг, (pi ) = 0 для <рч. G Нп. Заметим, что базисная

о

функция (pi , ассоциированная с вершиной X{j, принадлежит НП1 если G из. Следовательно, если все вершины треугольника лежат на dG, то для этого треугольника условие (23) выполнено. Для

о

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

Замечание 2. Несложно показать, что если на Тг коэффициенты kap постоянны и q = 0, то индикаторы коррекции 1\ и Ц точны.

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

При построении графа новой триангуляции будем исходить из требования сохранения конформности сетки. Второе важное требование: в процессе измельчения сетки не должно возникать треугольников с очень большими или малыми углами. Известно [22, 23], что появление треугольника с большим углом приводит к снижению точности конечно-элементной аппроксимации, а с малым — к росту числа обусловленности матрицы системы уравнений (5). Точность конечно-элементного решения зависит также и от гладкости сетки (соседние треугольники не должны резко отличаться друг от друга по площади).

Ниже описываются две процедуры измельчения сетки Ref-y и Ref2, отличающиеся реализацией первого шага этапа добавления новых узлов. Первая из них близка к процедуре, описанной в [24], а вторая является комбинацией метода бисекции [25] и указанной выше процедуры. В описании процедур используется понятие большей стороны треугольника, т.е. стороны, имеющей наибольшую длину. Если в треугольнике две или три стороны имеют максимальную длину, то в качестве большей стороны выбирается и фиксируется одна из них. I. Этап добавления новых узлов

1. Первый шаг. Обработка элементов списков Lists и Listt.

(a) Новые узлы ставятся на середины всех сегментов из списка Lists. Номера треугольников, соответствующих этим сегментам, заносятся во вспомогательный список List.

(b) Для каждого треугольника из списка Listt новые узлы ставятся или на середины всех сторон треугольника (процедура Refy), или на середину большей стороны (процедура Ref2). В первом случае в список List добавляются номера треугольников, смежных с треугольниками из Listt, во втором — номера всех треугольников из списка Listt.

2. Второй шаг. Для каждого треугольника из списка List выполняется следующая рекурсивная процедура.

(а) Если на середине большей стороны обрабатываемого треугольника нового узла еще нет, то он должен быть поставлен.

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

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

II. Этап построения графа новой триангуляции

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

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

Приведенное выше правило перестройки триангуляции позволяет избежать появления треугольников со слишком малыми или слишком большими углами.

5. Итерационный метод. Пусть и}^ — сетка на к-й итерации, а и^ — соответствующее приближенное решение задачи (1), (2) на этой сетке. Итерационный метод описывается следующим образом.

1. Инициализация. Задать начальную триангуляцию области С и сформировать пустой список Ыв^.

2. Построение начальной сетки. Построение сетки о?'0' осуществляется итерационно.

(a) Для каждого нового сегмента (вначале все сегменты считаются новыми) по формуле (7) вычислить значение К1 ив соответствии с описанным в пункте 3 правилом сформировать список ЫвЬд.

(b) Если список окажется пустым, то перейти к пункту 3, иначе, используя описанную в пункте 4 процедуру измельчения сетки, построить новую сетку и перейти к пункту 2а.

3. Начальный шаг. Построить систему уравнений (5), найти и'0' и по формуле (6) вычислить значение функционала /(и'0').

4. Итерационный шаг. Для к = 1,2,... выполнить следующие этапы.

(a) Для каждого нового треугольника Тг (вначале все треугольники считаются новыми) вычислить значение индикатора коррекции Р выбранного типа и по описанному в пункте 3 правилу отбора треугольников сформировать список Ыв^.

(b) Используя выбранную процедуру адаптивного измельчения сетки (см. пункт 4), построить новую сетку Построить систему уравнений (5), найти приближенное решение и^ и по формуле (6) вычислить ./(и^').

(c) Тест на сходимость. В случае выполнения условия

-J(uW) íC atol + rtol J(UW)

где atol и rtol — заданные допуски, итерационный процесс закончить, иначе к увеличить на 1 и перейти к пункту 4а.

Замечание. В тесте на сходимость можно использовать и другие критерии окончания итерационного процесса, в частности основанные на использовании апостериорных оценок погрешности zW = и^ — и*. Различные способы получения таких оценок рассмотрены в работах [12-15].

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

В качестве тестовой рассматривалась задача (1), (2), где С = {О ^ х\ ^ 1, 0 ^ х2 ^ 1},

2

кгг(х) = к22(х) = £- [е"*<* + е^"1'], г = 10"4,

а = 1

к12(х) = к21(х) = 0, д(ж) = е^и-^+М!-^

а /(ж) и иъ{х) были выбраны так, чтобы функция

2

и*(х) = 3+5] _ _ е(ха-1)/^"

а = 1

была классическим решением дифференциальной задачи (3), имеющим сильный экспоненциальный рост в приграничной области и в окрестности точки х = (7,7)- Для тестовой задачи с точностью Ю-13 имеем «/(и*) = -7,1573091002796.

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

В критерии отбора сегментов границы выбиралось ¡1 = 1/16 (в силу квадратичной зависимости погрешности интерполяции от шага при кусочно-линейной интерполяции гладкой функции на равномерной сетке) и полагалось е2и = Ю-5. В критерии отбора треугольников Тг использовалось значение в = 0,2.

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

Таблица 1 Метод на адаптивно измельчаемых сетках

к п т »V Г и) Гп

0 1424 573 280 4,58 • 10" 3 2,17 10" 1 2,37 • ю- 1

1 1616 661 296 1,67 • 10" 3 1,19 10" 1 1,86 • ю- 1

2 1752 721 312 9,97 • 10" 4 1,03 10" 1 1,19 • 10- 1

3 1960 805 352 5,02 • 10" 4 6,79 10" 2 1,16 • 10- 1

4 2424 1013 400 2,55 • 10" 4 3,19 10" 2 7,22 • 10- 2

5 2856 1181 496 1,60 • 10" 4 5,01 10" 2 5,59 • 10- 2

6 3632 1521 592 8,20 • 10" 5 3,19 10" 2 3,19 • 10- 2

7 3976 1693 592 6,88 • 10" 5 3,16 10" 2 3,16 • 10- 2

8 4824 2117 592 6,02 • 10" 5 1,47 10" 2 2,20 • 10- 2

9 6720 3065 592 4,99 • 10" 5 1,01 10" 2 2,19 • 10- 2

10 9752 4349 1056 3,42 • 10" 5 6,44 10" 3 7,85 • 10- 3

Таблица 2 Метод на равномерно измельчаемых сетках

к п т »V Г и) гп

0 18 4 12 5,78 • 10" 2 3,04 ю- 1 5,94 • ю- 1

1 72 25 24 1,61 • 10" 2 2,09 10- 1 5,46 • ю- 1

2 288 121 48 5,65 • 10" 3 1,71 10- 1 4,67 • ю- 1

3 1152 529 96 1,71 • 10" 3 1,19 10- 1 3,46 • ю- 1

4 4608 2209 192 4,03 • 10" 4 5,36 10- 2 2,09 • ю- 1

5 18432 9015 384 9,51 • 10" 5 1,58 10- 2 9,74 • ю- 2

6 73728 36481 768 2,33 • 10" 5 3,62 10- 3 3,73 • ю- 2

7 294912 146689 1536 7,87 • 10" 6 8,98 10- 4 1,21 • ю- 2

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

Рис. 1. Триангуляция области (а) и начальная сетка (б)

Рис. 2. Сетка, построенная на последней итерации (а), и полученное на ней решение (б)

В таблицах (индекс к опущен) использованы следующие обозначения:

и^(х)-и*(х) ь и^(х)-и*(х) ь ,7{и^)-,1(и*)

гш = тах

и*(х)

к

= тах

хвО.(к)

и*(х)

Г$ =

|7И|

Лк)

где ПМ = у п\к), а Щ

(к)

множество узлов вспомогательной сетки в треугольнике Тг. Использо-

¿=1

вание множества позволяет оценить погрешность в большом числе точек, расположенных между узлами сетки и тем самым смоделировать оценку для погрешности в норме пространства С (Сг). Задача считалась решенной, когда относительная погрешность полученного решения г^ не превышала 10"2.

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

Эксперименты показали, что наиболее эффективным для рассматриваемой тестовой задачи оказался метод, использующий адаптивную стратегию S2 и процедуру Ref2 измельчения сетки. Стандартный метод конечных элементов с равномерным измельчением сетки существенно уступает любому из вариантов (комбинации стратегий и процедур измельчения сетки) предложенного в работе метода как по требуемой памяти, так и по времени решения задачи. Он требует примерно в 25 раз больше узлов и в 10 раз больше времени, чем наилучший адаптивный метод, предложенный в данной работе.

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

1. Babuska I.,Rheinboldt W. С. A posteriory error estimates for the finite element method / / Intern. J. Numer. Meth. Engrg. 1978. 12. N 11. P. 1597-1615.

2. Babuska I., Rheinboldt W. C. Error estimates for adaptive finite element computations / / SIAM J. Numer. Anal. 1978. 15. N 4. P. 736-754.

3. Bank R. E., Weiser A. Some a posteriory error estimates for elliptic partial differential equations // Math. Comput. 1985. 44. N 170. P. 283-301.

4. Ainsworth M., О den J.T. A procedure for a posteriory error estimation for h-p finite element methods // Comput. Meth. Appl. Mech. Engrg. 1992. 101. N 1. P. 73-96.

5. Ainsworth M.,Oden J. T. A unified approach to a posteriory error estimation using element residual methods // Numer. Math. 1993. 65. N 1. P. 23-50.

6. Bornemann F.A., Erdmann B.,Kornhuber R. A posteriory error estimates for elliptic problem in two and three space dimensions // SIAM J. Numer. Anal. 1996. 33. N 3. P. 1188-1204.

7. Zienkiewicz O.C., Kelly D.W., de Gago J.P., Babuska I. Hierarchical finite element approaches, adaptive refinement and error estimates // The Mathematics of Finite Elements and Applications / Ed. J.R. Whiteman. Academic Press, 1982. P. 313-346.

8. Zienkiewicz O.C., Craig A. Adaptive refinement, error estimates, multigrid solution, and hierarchical finite element method concepts // Accuracy Estimates and Adaptive Refinements in Finite Element Computations / Eds. I. Babuska, О. C. Zienkiewicz, J. P. de Gago, E. R. de Oliveira. John Wiley & Sons, 1986. P. 25-59.

9. Bank R.E.,Smith R.K. A posteriory error estimates based on hierarchical bases // SIAM J. Numer. Anal. 1993. 30. N 4. P. 921-935.

10. Николаев E.C. Метод решения краевой задачи для ОДУ второго порядка на последовательности адаптивно измельчаемых и укрупняемых сеток // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2004. № 4. С. 5-16.

11. Николаев Е.С. Метод решения краевых задач для систем линейных ОДУ на адаптивно измельчаемых сетках // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2000. № 3. С. 11-19.

12. Verfiirth R. A posteriori error estimation and adaptive mesh-refinement techniques //J. Comput. Appl. Math. 1996. 50. N 1-3. P. 67-83.

13. Ainsworth M., Oden J.T. A posteriory error estimation in finite element analysis // Comput. Meth. Appl. Mech. Engrg. 1997. 142. N 1. P. 1-88.

14. Jones M.T., Plassman P. E. Adaptive refinement of unstuctured finite-element meshes // Finite Elem. Anal. Design. 1997. 25. N 1-2. P. 41-60.

15. Zhu J.Z. A posteriori error estimation — the relationship between different procedures // Comput. Meth. Appl. Mech. Engrg. 1997. 150. N 1-4. P. 411-422.

16. Huang W., Ren Y., Russel R. D. Moving mesh partial differential equations (MMPDES) based on the equidistribution principle // SIAM J. Numer. Anal. 1994. 31. N 3. P. 709-730.

17. White A.B. On selection of equidistributing meshes for two-point boundary-value problems // SIAM J. Numer. Anal. 1979. 16. N 3. P. 472-502.

18. Дегтярев Л.М., Дроздов В.В., Иванова Т.П. Метод адаптивных к решению сеток в сингулярно-возмущенных одномерных краевых задачах // Диф. ур-ния. 1987. 23. № 7. С. 1161— 1168.

19. Felippa С. A. Optimization of finite element grids by direct energy search // Appl. Math. Model. 1976. 1. N 1. P. 93-96.

20. Delfour M.,Payre G.,Zolesio J.-P. An optimal triangulation for second-order elliptic problems// Comput. Meth. Appl. Mech. Engrg. 1985. 50. N 3. P. 231-261.

21. To u rig n у Y., Hiilsemann F. A new moving mesh algorithm for the finite element solution of variational problems // SIAM J. Numer. Anal. 1998. 35. N 4. P. 1416-1438.

22. Babuska I., Aziz A. K. On the angle condition in the finite element method // SIAM J. Numer. Anal. 1976. 13. N 2. P. 214-226.

23. Fried I. Condition of finite element matrices generated from nonuniform meshes // AIAA J. 1972. 10. N 2. P. 219-221.

24. Verfiirth R. A posteriori error estimators for the Stokes equation non-conforming discretization // Numer. Math. 1991. 60. N 2. P. 235-249.

25. Rivara M.-C. Mesh refinement processes based on the generalized bisection of simplices // SIAM J. Numer. Anal. 1984. 21. N 3. P. 604-613.

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

УДК 519.6

M. Дана, Х. Д. Икрамов

О РЕШЕНИИ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ, МАТРИЦЫ КОТОРЫХ ЯВЛЯЮТСЯ МАЛОРАНГОВЫМИ ВОЗМУЩЕНИЯМИ ЭРМИТОВЫХ МАТРИЦ

(кафедра общей математики факультета ВМиК)

1. Введение. Пусть нужно решить систему линейных уравнений

Ах = Ь (1)

с нормальной га X га-матрицей А. Нормальность в данном контексте означает выполнение условия

АА* = А* А. (2)

Для разреженной неэрмитовой матрицы А стандартным выбором алгоритма для решения системы (1) является обобщенный метод минимальных невязок GMRES (см. [1, раздел 6.6.6]). Однако в этом методе не удается использовать специфику матрицы А, выражаемую равенством (2): алгоритм Арнольди, лежащий в основе GMRES, приводит нормальную матрицу А к той же форме Хессенберга, что и матрицу общего вида. В результате сохраняются два главных недостатка, присущие GMRES по сравнению с методами типа Ланцоша: большой расход памяти и зависимость количества арифметической работы, выполняемой на отдельном шаге, от его номера к.

При дополнительном предположении, что спектр матрицы А принадлежит алгебраической кривой Г невысокого порядка (примером реалистической ситуации, где это предположение выполнено, могут служить нормальные матрицы с почти-вещественным спектром), в [2] предложен иной подход к решению системы (1), который вкратце состоит в следующем. Методы типа Ланцоша и GMRES ищут приближенные решения в крыловских подпространствах возрастающей размерности, т.е. в линейных оболочках конечных отрезков степенной последовательности

х, Вх, В2х,... . (3)

Здесь В — это или сама матрица А, или некоторая матрица, тесно с А связанная; в качестве х обычно (при отсутствии хорошего начального приближения) берут правую часть Ь системы (1) или полученной из нее "предобусловленной" системы.