Применение энергетических методов при построении векторных полей Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 9, № 4, 2004

ПРИМЕНЕНИЕ ЭНЕРГЕТИЧЕСКИХ МЕТОДОВ ПРИ ПОСТРОЕНИИ ВЕКТОРНЫХ ПОЛЕЙ*

В. В. ДЕНИСЕНКО Институт вычислительного моделирования СО РАН,

Красноярск, Россия e-mail: denisen@icm.krasn.ru http://icm.krasn.ru/personal.php?persid=38

A principle of minimum of a quadratic energy functional is formulated for 3D problem of vector field calculation from given curl and divergence in the domain with impenetrable boundary. The obtained matrix of finite element method is shown to be symmetric. It is proved to be positively defined for a number of particular cases.

Введение

Построение векторного поля по заданным плотностям вихря G и источника Q требуется при решении многих задач механики и физики [1]. Проще строятся безвихревые поля, например электростатическое поле в вакууме E или поле скорости при безвихревом течении, для которых G = 0. В таких случаях принято вводить новую неизвестную функцию — потенциал E = -grad при нахождении которого приходится решать те или иные краевые задачи для уравнения Пуассона. Эффективным методом построения приближенных решений задач для уравнения Пуассона является энергетический [2]. К такой же задаче с иными правыми частями исходная задача может быть сведена и при ненулевом роторе, если построить некоторую вспомогательную векторную функцию с заданным ротором. Универсальным способом построения такой функции является интегрирование закона Био — Савара, но он требует объема вычислений, на порядки превосходящего затраты на решение задачи для <^.

Для общего случая, когда среда является неоднородной, анизотропной и трансверсально анизотропной, Q = 0, G = 0, энергетический метод для трехмерных задач был предложен и обоснован в [4], для двумерных — в [3].

Для полей с нулевой дивергенцией, например для магнитного поля B, принято вводить векторный потенциал A , такой что B = rot A. Задание нормальной компоненты B на границе означает задание ротора касательных компонент A. Энергетический метод можно применить, если предварительно преобразовать это граничное условие, решив вспомогательную двумерную задачу на граничной поверхности. Следует построить векторную функцию, имеющую только касательные к границе компоненты, с заданным ротором и

* Работа выполнена при поддержке Российского фонда фундаментальных исследований (грант № 0405-64088) и программы ОФН-16 РАН (проект 2.15.2).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

нулевой дивергенцией. Необходима их продолжимость в область с конечной энергией аналогично тому, как это необходимо для неоднородной задачи Дирихле для уравнения Лапласа [5]. Тогда граничное условие заменяется на задание касательных компонент А.

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

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

1. Формулировка задачи

В стационарном случае при единичной магнитной проницаемости магнитная индукция B = (Bx , By ,Bz) удовлетворяет уравнениям

rot B = G, (1)

div B = Q, (2)

где в случае магнитостатики заданная функция Q = 0, но мы используем общий вид правой части, поскольку в задачах для других векторных полей функция Q отлична от нуля, например, в задачах электростатики это заданная плотность заряда. Здесь x,y,z — декартовы координаты, плотность стороннего тока G — заданная функция координат, которая для разрешимости задачи должна иметь нулевую дивергенцию

div G = 0, (3)

поскольку дивергенция левой части уравнения (1) тождественно равна нулю.

Предполагаем, что область П — ограниченная и ее граница Г — гомеоморфная сфере дважды непрерывно дифференцируемая поверхность. Нормальную и касательные к границе компоненты векторов будем отмечать индексами n и l.

Если за границей рассматриваемой области находится идеальный проводник, то на

границе нормальная компонента магнитной индукции B обращается в нуль. Такое же

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

Bnlr = q, (4)

(q — заданная функция).

Для разрешимости задачи заданные функции Q и q должны удовлетворять условию

J Q dn = j) qdr. (5)

В противном случае получим противоречие, проинтегрировав равенство (2) по всей области и (4) — по всей границе, поскольку интегралы левых частей равны в силу теоремы Гаусса — Остроградского.

При сделанных предположениях решение задачи (1), (2), (4) единственно. Например, в [4] это доказано из энергетических соображений.

2. Краевая задача для пары потенциалов

Введем скалярный F и векторный P потенциалы так, что

B = -grad F + rot P, и потенциалы удовлетворяют условиям

Pilr = 0; (6)

J FdQ = 0. (7)

Здесь Pi — обе касательные к границе компоненты. Как следствие,

rotra P|r = 0.

Для устранения неоднозначности P используем условие

div P = 0. (8)

Из исходных уравнений (1), (2), (4) получаем уравнения для введенных потенциалов:

rot rot P = G; (9)

-div grad F = Q; 3F_ dn

= q. (10)

г

Задачи для F и P расщепились: (7), (10) — для F и (6), (8), (9) — для P.

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

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

Поэтому ниже будем рассматривать только задачу (6), (8), (9) для P. Соответствующую P часть B отметим индексом P. Она дает полное решение исходной задачи, если Q = 0, q = 0, поскольку в этом случае задача (7), (10) имеет только тривиальное решение F = 0.

Решение задачи (6), (8), (9) можно свести к минимизации функционала энергии

W (P) = 2[P, P] -J PT G dn (11)

на множестве функций, удовлетворяющих условиям (6). Квадратичная часть W(P) соответствует энергетическому скалярному произведению

[v, P] = / ((rot v)T rot P + div v div P) dn. (12)

Эта симметричная билинейная форма на рассматриваемом множестве функций положительно определена. Энергетическая норма эквивалентна сумме норм декартовых компонент Р как элементов пространства Ж2(1)(П). Минимум функционала энергии существует, единствен и дает обобщенное решение задачи (6), (8), (9), которое при дополнительном предположении гладкости является классическим решением. Доказательства для более общего случая проведены в [4].

С учетом аналогичных свойств задачи для Г принцип минимума функционала энергии справедлив и для всей исходной краевой задачи (1), (2), (4).

Квадратичная форма, соответствующая энергетическому скалярному произведению

(12), может быть интегрированием по частям преобразована к виду

[Р, Р] = У Рк)2 СП + У Р + [Р х гоі Р]п - £ Ри

СГ, (13)

где к — номер декартовой компоненты.

Поскольку рассматриваются только функции Р, удовлетворяющие (6), векторное произведение в интеграле по границе имеет нулевую нормальную компоненту. С учетом этого же условия (6) оставшаяся разность равна

Р 1д. (БР ) _ р дРП = Р21 дБ (14)

Рп Бди( п) п дп = Рп Б дп, (14)

где Б — площадь параллельной границе малой площадки, ограниченной нормалями к границе. На плоских участках границы Б = 1, на цилиндрических — Б = г, на сферических — Б = г2. В общем случае Б есть произведение соответствующих коэффициентов Ламэ.

После проведенного преобразования квадратичной формы функционал энергии (11) принимает вид

И'"(Р) = £ / (|(«га^ Рк)2 - РкСк) СП + 2 ^ Р2БдПсг. (15)

Заметим, что интеграл по области расщепился на три интеграла для трех декарто-

вых компонент векторной функции Р, причем каждый интеграл имеет тот же вид, что и для скалярного потенциала. Отличие от независимой минимизации трех функционалов определяется наличием интеграла по границе в (15), куда входят все три компоненты Р.

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

3. Численная реализация граничного условия

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

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

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

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

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

После покомпонентных вычислений получаем интеграл по треугольнику вида

При этом не только Рх, но и компоненты Рх, Ру отличны от нуля.

Поскольку производные по х, у кусочно-линейных функций постоянны в треугольнике, их можно вынести за знак интеграла, а интегралы от самих кусочно-линейных функций есть средние арифметические их узловых значений, умноженные на площадь треугольника. Пронумеруем вершины треугольника в направлении против часовой стрелки і = 1, 2, 3 и обозначим через г их координаты, через п — определенные в них нормали к границе, через РП — значения нормальной компоненты функции Р в этих узлах. Тогда значения функции Р в этих узлах

Это условие в совокупности с линейной интерполяцией декартовых компонент Р является способом аппроксимации главного краевого условия (6).

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

Все векторные произведения обращаются в нуль, если все три нормали параллельны.

Сумма квадратичных форм вида (18), соответствующих всем граничным треугольникам, должна быть добавлена к квадратичной и линейной формам, получающимся из интеграла по объему в (15). Последние определены для произвольных значений компонент Р в граничных узлах, но нас они интересуют только для Р^ вида (17).

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

(16)

(17)

(18)

Зафиксируем г-й узел, лежащий на границе области, и просуммируем квадратичные формы (18) по прилегающим к нему треугольникам. Вершины ломаной, ограничивающей эту совокупность треугольников, пронумеруем индексом ] = 0, 1, ^0, где 0 и к соответ-

ствуют одному и тому же узлу, поскольку ломаная замкнута. Получаем

1 30

- рпх; ["<х п з - г,-1>^.

3 = 1

Производная этой квадратичной формы по РП равна

1 30

- Дп X Пз](Гз+1 - Г,-_1)РП. (19)

3=1

Обозначим через О,, О;, О; производные от интеграла по объему в (15) по Р,, Р;,Р; соответственно. Здесь г — номер произвольного граничного узла сетки. Поскольку свободным параметром в этом узле является только Ргп ив силу (17)

РХ = 4р\ р; =п; рп, р; = п; рп, (20)

для граничного узла получается одно (а не три) уравнение, соответствующее равенству

нулю производной функционала энергии (15) по Рхп. Используем правило дифференцирования сложных функций и добавим вклад от интеграла по границе (19) с соответствующим (15) множителем 1/2. Получаем уравнение

дШ ........................1 ^

дРп = пХ ох + п; о; + п; о; + - £[щ х п, ], - г^орп = о. (21)

х 3=1

В каждом таком уравнении ОгХ есть линейная комбинация значения РХ, значений Р,, соответствующих соседним граничным узлам ] = 1,]0, а также значений РХ в тех внутренних узлах, которые соединены одним ребром сетки с рассматриваемым граничным узлом г. При этом следует все граничные значения Р’Х, Р, выразить через Рхп, Рп по формулам (20). Аналогично преобразовав О;, О;, получаем из (21) равенство нулю линейной комбинации РХП, РП плюс определенной линейной комбинации значений Р в приграничных внутренних узлах и плюс член пт ОХ, который получится после умножения свободных членов Юх на пх.

Если зафиксировать значения Р во внутренних узлах, система уравнений (21) представляет собой систему г0 уравнений для г0 неизвестных Рхп . Матрица этой системы симметрична по построению. Поскольку аппроксимирующие кусочно-линейные функции не удовлетворяют поточечно-главным краевым условиям (6), из положительной определенности квадратичной части функционала энергии, доказанной в энергетическом пространстве, не следует положительная определенность этой матрицы.

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

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

4. Положительная определенность матрицы вариационно-разностной схемы

Пусть область П — прямой круговой цилиндр единичного радиуса и высотой г0. Построим на боковой поверхности равномерную сетку ах = гка, г, = ^'к^, где г = 1,, ] = 0,3. Здесь г, а, г — цилиндрические координаты; положительные числа ка, к;, кг — шаги сетки. Полагаем I четным и

где число к — шаг сетки и константы не зависят от к. Удобно иметь сетку с регулярной нумерацией узлов, такой же, как в параллелепипеде. Для использования кусочно-линейных функций следует каждую кубическую в смысле нумерации ее узлов ячейку сетки разбить на тетраэдры. Мы делаем это, определив некоторый центр ячейки и центры ее шести граней. Тогда каждая грань естественным образом делится на четыре треугольника, а ячейка — на 24 тетраэдра. Значения декартовых компонент Р в этих центрах получаем некоторой интерполяцией узловых значений, за исключением граней, аппроксимирующих границу. Каждый такой сеточный прямоугольник разделим диагоналями на четыре треугольника. Построим в центре перпендикуляр высотой кг, который будет общим ребром для четырех тетраэдров, прилегающих к граничным треугольникам. Считаем, что разбиение области на тетраэдры выполнено с соблюдением стандартных ограничений для нерегулярной сетки, которые обеспечивают эквивалентность Ь2 нормы кусочно-линейной функции и евклидовой нормы набора ее узловых значений:

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

В центре каждой прямоугольной ячейки на боковой поверхности, который нумеруем полуцелым индексом і + 1/2, і + 1/2, нормаль пі+1/2а-+1/2 строим усреднением из четырех равноценных вершин, поэтому пі+1/2;^+1/2 = (1, 0, 0).

В силу (17) значение Р в граничном узле определяем одним параметром Рп. Значение Рн-і/2 і+1/2 во вспомогательном узле строим усреднением из четырех узлов

ка = Са к, кг = Сг Л,, кг = сг Л,,

(22)

(23)

(24)

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

Поскольку узловые значения Р- обращаются в нуль при ] = 0 и ] = 3 и периодичны по г, их можно представить в виде суммы:

Pn

Pij

2

J _1

sin (lzjn/z0)

і/2

^(aim sin (mai) I bim cos (mai)) I bro/v^

m=1

(25)

Считаем aij/2 = 0, b0,i/2 = 0, поскольку это коэффициенты перед тождественно равной нулю функцией sin (in), и ai,0 = 0, поскольку это коэффициент перед sin0.

Несложно проверить, что все функции, по которым проведено разложение (24), взаимно ортогональны и имеют единичные нормы, если в качестве скалярного произведения использовать

I J

(P^) = ^^ p” р

Имеем также соотношение

I J

ijij

i=1 j=0

J _1 I/2

££(PJ)2 = ££(a?m I bm) •

i=1 j=0

(26)

i=1 m=0

Для любого треугольника, лежащего на торцах цилиндра, выражение (18) обращается в нуль, поскольку или все три нормали параллельны, или один из узлов лежит на окружности г = 1, но в этом случае Рп = 0. Поэтому в граничном интеграле из (13) можно оставить только треугольники из Г^. Для четверки треугольников, образующих прямоугольник, используя явные выражения нормалей и координат узлов, можно записать (18) в виде

hz ЙІП К Рп Рп , К ЙІП К Рп Рп , 4^ 8Іп(^«/2)(рп ,2 (27)

6 Рі,і Рі+1^ + 6 Рг,і+1Рг+1,і+1 + 3 (Рг+1/2,і+1/2) • (27)

Pi

ij

Просуммируем по всем прямоугольникам, подставив выражение (24). Поскольку 0 при ] = 0 и при ] = 3 и есть периодичность по г, суммы по г,некоторых

слагаемых из (27) равны. Получаем

h

-z (sin ha I sin (ha/2))Y^ Pnj Pi+1,jI

і,]

h- sin (ha/2)

3

n 2 n n n n n n

/ , l(Pi,j ) + pi,j Pi,j+1 ^ pi,j Pi+1,j + 1/2 ^ Pi,j + 1Pi+1,j /2J •

і,]

Выполнив суммирование по г,^ с учетом ортонормированности функций, получаем следующее выражение для интеграла по границе из (13):

hz sin (ha/2)

З

J_1 I/2

££

i=1 m=0

2 cos -2 I 1^ cos (mha) I 1 I cos----1

I cos (mha) cos

lnhz

Z0

(28)

I

Поскольку тЛ,а изменяется от 2п/1 до п при т = 1,...,//2, эта квадратичная форма не является положительно определенной. Отрицательны слагаемые, соответствующие мелкомасштабным гармоникам. Независимо от I при т > 1/3 выражение в квадратичных скобках отрицательно, а при т < 1/4 положительно.

Наряду с рассмотренной суммой по всем граничным треугольникам выражений (18), которая является аппроксимацией интеграла по границе в (13), в квадратичную форму

(13) еще входит интеграл по области.

С помощью легко проверяемого неравенства

(grad |P|)2 < 3 ^(grad Pk)2

(29)

k=1

можно оценить интеграл по области снизу.

В силу (17) в граничных узлах |Pj | = pj. Внутри каждого приграничного тетраэдра модуль вектора grad не меньше модуля его касательной к граничному треугольнику составляющей. Поэтому для четверки тетраэдров

(grad |P|)2dQ >

ha hz hr

24

pn Pn\ 2 / pn pn \ 2

pi+1,j pij \ + / pi+1,j+1 pi,j+1 \ +

ha

hn

P»n ]Dn\ 2

p

' pn pn \ 2

+ I \ + I pi+1,j+1 pi+1,j

hz

hz

Просуммировав это неравенство по всем приграничным тетраэдрам, усилим его, распространив интегрирование неотрицательной функции на все тетраэдры, и учтем периодичность pj по i и pj = 0 при j = 0 и j = J. Получаем

/(grad |р|)2Ж > £(pn+1 j - pij)2 + £( j - pi)2.

nh a ij z ij

Подставим в правую часть выражения pj вида (25) и просуммируем по i, j с учетом ортонормированности функций. С учетом (29) получаем оценку снизу для интеграла по области из (13):

^ I (grad pk )2 dQ fc=1Oh

I 1/2

2 JO ^ hahzhr ((sin2 (mha/2) (sin2 (/nhz/(2zo)U 2 2

k) ^--- 2^ 2^ 1 --------h2-------- + ------------------------ ) (a1m + blm)-

3

l=1 m=0

hz2

Прибавив выражение (28) граничного интеграла из (13), получаем для квадратичной формы оценку снизу

sin2(mha/2) sin2(/nhz/rz0)

h2

+

hz2

sin(ha/2) /Y0 ha 1nh^\ ^nh2

+----iTT---- I I 2 cos "7Г + 1 + cos---- cos(mha) + 1 + cos I -------

hahr \ \ 2 z0 / \ z0

(alm + bim) •

3

f

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

[P P] г ha hz hr

з

2sin(ha/2) ha hr

^ _ 4sin(ha/2)(1 + cos(ha/2))

ha hahr

+ l к-

1 4sin(ha/2),sm2((nhz/(2zo))

hahr

(cos(ha/2) + 2) +

sin2(mha/2) +

(aim + bim) •

Еще усилим неравенство, используя 0 < cos (ha/2) < 1 и ha/4 < sin (ha/2) < ha/2:

[P, P] г

hahz h,

aibz1 EE

І m

з

hr + 1 ha /ir *sin

тЛ0

_2_

+ ' й- hr|sin

/nhz

2zn

(aim+bim) •

Поскольку шаги сетки по разным направлениям должны быть одного порядка малости (22), можно считать, что Лг > Ла(4Ла),ЛГ > (2^). Поэтому коэффициенты перед

квадратами синусов неотрицательны и можно усилить оценку:

[p. p] г

С учетом (26) и (22) получаем

£(«L + О-

[p. p] г ^ h^ (P" )2

(ЗО)

i,,

Поскольку коэффициенты матрицы системы алгебраических уравнений (21) для Рг” есть вторые производные квадратичной формы (13) по Р", Р™^,, полученная оценка (30) означает положительную определенность матрицы системы (21).

Теперь докажем положительную определенность всей матрицы вариационно-разностной схемы.

Обозначим через У совокупность всех значений Рг”, через X — всех значений трех декартовых компонент Р во внутренних узлах. Тогда квадратичную форму (13) можно записать как квадратичную форму с блочной структурой матрицы:

[Р, Р] = (Xт ,Ут) Неравенство (30) означает, что

т ^TW A B

BT C

X

Y

XTAX + 2XTBY + YTCY г h2YTY.

г 3

(З1)

(З2)

В подпространстве У = 0, т. е. на функциях Р, все три компоненты которых в граничных узлах равны нулю, положительная определенность следует из положительной определенности (12), поскольку условию (6) такие функции удовлетворяют поточечно. Имеем

XT AX >

1

С2

p2 da

1

2

и в силу (23)

ХтАХ > — Л3ХТX. (33)

С2

Из (32) следует положительная определенность матрицы С, которая и является матрицей системы (21):

утсу > 2£а££ л2УтУ. (34)

- 3 V у

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

В = v/A"тBv/C"т. (35)

Как и всякую квадратную матрицу, В можно представить в виде произведения ортогональной матрицы V на симметричную неотрицательно определенную матрицу К:

В = уК, (36)

где К находится из условия

К2 = Вт В. (37)

Минимум квадратичной формы из (31) по X достигается при таком X*, что

АХ* + ВУ = 0,

т. е. при

X* = -А-1 ВУ.

Вычислив соответствующее такому X* минимальное значение левой части в (32), получаем неравенство

Ут(С - ВтА-1В)У > 2СаС^Л2УтУ,

справедливое при любом У и поэтому означающее некоторую ограниченность матрицы В. Перепишем левую часть, используя (35), (37):

Ут^С(I - К2)УСУ > 2С^Л2УтУ.

3

Здесь I — единичная матрица.

Обозначим ^ = -\/СУ:

^т(I - К2)^ > 2Са^Л2^тС"1^. (38)

3

Чтобы продолжить оценку снизу, рассмотрим матрицу С. Элемент матрицы Сгт получается как энергетическое скалярное произведение двух кусочно-линейных векторных функций Р;, Рт, равных нулю во всех узлах, кроме Р™ = 1 для одной функции и Рт = 1 для другой.

Используя для квадратичной формы (31) исходное выражение (12), получаем

С;

1т

J(&уР; &уРт + гс^Рг го1Рт)^П.

В силу неравенства Коши — Буняковского

С2т < ^ ((&уРг)2 + (хОРг)2)dQ J((&уРт)2 + (го1Рт)2^П

Выберем в качестве I тот узел, для которого такой интеграл максимален. Тогда

|Си</ ((^уРг )2 + МРг )2^П.

Индекс I при оценке этого интеграла опускаем.

Воспользовавшись легко проверяемым неравенством

получаем

2

dП.

(39)

В соответствии с приведенными в начале раздела правилами разбиения ячеек сетки на тетраэдры значение Р™ = 1 при нулевых значениях во всех остальных узлах интерполируется в центры восьми четырехугольников, |Р| = 1/4, и в центры четырех ячеек, |Р| = 1/8. Получаются 96 тетраэдров, в которых кусочно-линейная функция Р не обращается в нуль тождественно. Они делятся на три группы. В вершинах 24 тетраэдров первой группы |Р| = (1,1/4,1/8, 0). В вершинах 24 тетраэдров второй группы |Р| = (1/4,1/8, 0,0). В вершинах 48 тетраэдров третьей группы |Р| = (1/8, 0, 0, 0). В тетраэдрах первой группы функцию Р можно рассматривать как сумму трех функций, каждая из которых не равна нулю только в одной из вершин. В силу неравенства Коши — Буняковского интеграл из (39) по тетраэдру не превосходит величины

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

где через с5Л, обозначена минимальная высота тетраэдров.

Для двух оставшихся функций, Р(2) и Р(3), аналогично получаем оценки с константами, в 16 и 64 раза меньшими. Поэтому для тетраэдров первой группы подынтегральное

выражение в (40) не превосходит 13/(12е§^2). Аналогично рассмотрев тетраэдры второй и третьей групп, получаем меньшие значения. Поскольку суммарный объем всех 96 тетраэдров не превосходит 4^^Пг, с учетом (39) и обозначений (22) получаем

I Г1 | ^ 39агг ,

|С1т,| < 2

с5

Если узел т не является соседним с узлом /, то Сгт = 0 и, значит, сумма |Ст| по всем т не более чем в 9 раз превосходит максимальное значение | Сгт |:

351^сасг сг/с;2. (41)

Собственные числа матрицы не превосходят суммы по строке модулей ее элементов. Поэтому максимальное собственное число матрицы С не превосходит значения (41), т. е.

уТСУ < 351с2^Сг ПУТУ (42)

С5

и, значит,

с2

С-1^ > ------5--- Z.

351сасг сг п

Подставляя эту оценку в неравенство (38), получаем

с2

ZT(I - К2^ > —ZTZ.

1 ; “ 528сг

Коэффициент в правой части положителен. Он не может быть больше 1. Это легко проверить, рассмотрев левую часть в базисе, где К становится диагональной матрицей. Поэтому его можно обозначить через с6(2 — с6), где 0 < с6 < 1:

ZT(I — К2^ > Сб(2 — се^тZ.

Неотрицательная определенность К2 и это неравенство для самой матрицы К означают для любого Z

ZT(I — К^ > C6ZT^ (43)

что тоже легко установить, приведя К к диагональному виду.

Воспользуемся неотрицательностью квадрата вектора:

(\ГКут\[ах + л/К/СУ)т(\ГКут\[ах + л/К/СУ) > 0,

откуда следует

2|хт\/А/к/\/Су| < хт\/а/Ку^уАх + /^К^у.

С учетом обозначений (35)-(37) получаем

2|хтву| < (УАх)туКVт(УАх) + (Усу)тК(^су).

Используем это неравенство для оценки квадратичной формы (31), записав квадратичные формы X и У как квадраты \J~AX и ^У:

хтАх + 2хтву + утсу > (УАх)т(/ — уК/т)(^Ах) + (Усу)т(/ — К)(УСу).

Поскольку У — ортогональная матрица, с учетом (43) получаем неравенство

х т ах + 2х т СУ + У т СУ > с6(х т Ах + У т СУ),

которое с учетом положительной определенности матриц А (33) и С (34) и положительности с6 означает положительную определенность всей матрицы вариационно-разностной схемы:

хтах + 2хтсу + утсу > с^ С2 п3хтх + П2УтУ^ .

Напомним, что константы не зависят от п, они определяются только ограничениями на форму сеточных тетраэдров, но не их размером.

При достаточно малых П, П < 2сасгс2/(3с3), эту оценку можно усилить:

(х т ,У т ) ( вАт С)(У) > с7п3(х т ,У т ) ( У). (44)

Аналогично (42) может быть получена верхняя оценка собственных чисел всей матрицы из (31) с иной константой, но тоже первого порядка по п. После этого с учетом оценки снизу (44) получаем ту же оценку 0(1/П2) для числа обусловленности матрицы вариационно-разностной схемы, что и для классической семиточечной разностной схемы уравнения Пуассона в параллелепипеде.

Заключение

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

При аппроксимации энергетического пространства кусочно-линейными функциями получается вариационно-разностная схема с симметричной матрицей.

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

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

Список литературы

[1] СЕДОВ Л.И. Механика сплошной среды. Т. 2. М.: Наука, 1970. 568 с.

[2] Михлин С.Г. Вариационные методы в математической физике. М.: Гостехиздат, 1957. 378 с.

[3] ДЕНИСЕНКО В.В. Энергетические методы для эллиптических уравнений с несимметричными коэффициентами. Новосибирск: Изд-во СО РАН, 1995. 204 с.

[4] ДЕНИСЕНКО В.В. Энергетический метод для трехмерных эллиптических уравнений с несимметричными тензорными коэффициентами // Сиб. мат. журн. 1997. Т. 38, № 6. С. 1267—1281.

[5] ГОДУНОВ С.К. Уравнения математической физики. М.: Наука, 1979. 392 с.

Поступила в редакцию 18 февраля 2004 г., в переработанном виде —16 апреля 2004 г.