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

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

Том 14, № 1, 2009

Оптимизация структуры матриц массы и жесткости для векторного конечно-элементного базиса на ортогональных сетках

М.В. Пузлнов, Э.П. Шурина

Новосибирский государственный технический университет, Россия e-mail: misha. puzanov@gmail. com

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

Ключевые слова: уравнения Максвелла, векторный метод конечных элементов, ортогонализация базиса.

Введение

Многие нестационарные задачи вычислительного электромагнетизма при использовании конечно-элементной дискретизации в пространстве H(rot, Q) приводят [1, 2] к решению системы обыкновенных дифференциальных уравнений вида

где u = aiNi; Nj £ H(rot, П) — векторный базис, используемый при дискретизации;

i= 1

матрицы R, M1, M2 — квадратные, размерности п. В случае, если u соответствует полю E, то

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

Использование неявных схем по времени (схемы Ньюмарка, Кранка—Николсона, Адамса—Молтона и др.) теоретически приводит к безусловно устойчивому вычислительному процессу, который требует решения СЛАУ (системы линейных алгебраических уравнений) с матрицей вида аК + ЪЫ£ + оЫ3, а,Ъ,е £ К, на каждом шаге по времени. В случае малых временных шагов для достижения заданной точности при

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

(1)

n

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

Известно, что использование явной схемы по времени приводит к условно устойчивому вычислительному процессу, в котором устойчивость определяется соотношением Куранта. Применение явной схемы требует решения СЛАУ с матрицей вида pM£ + qMs, p,q Е R, на каждом шаге. Таким образом, уменьшение меры обусловленности матрицы массы M, Mij = J Nj • Nj и количества ненулевых элементов в ней способствует

п

снижению затрат времени вычислений и необходимой памяти.

Векторный МКЭ [3 - 6] предоставляет некоторые возможности конструирования базисов (в H(rot, П) и H(div, П)), приводящих к уменьшению числа ненулевых элементов в глобальных матрицах жесткости и массы и обеспечивающих эффективное обращение матрицы массы.

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

В [7, 8] предлагается метод приведения матрицы массы к диагональному виду заменой точного интегрирования при вычислении скалярного произведения приближенной процедурой, использующей квадратуры Гаусса—Лежандра или Гаусса—Лобатто. Все полиномы, формирующие базисные функции, ортогональны на системе узлов интегрирования. Такая процедура приводит к некоторой потере точности, которая, однако, может быть нивелирована достаточно высоким порядком полиномов. Основное преимущество данного метода в том, что он применим и к неортогональным гексаэдральным сеткам.

Другой способ уменьшения затрат на обращение матрицы массы — ее приближенное обращение [9], в котором из матрицы K = M-1 удаляются элементы, удовлетворяющие

критерию |fcj7-1 < -s——. Этот вариант, соответствующий областям с однородными

minq Kqq

физическими свойствами, может давать плохую аппроксимацию M-1 в случае сильно меняющихся коэффициентов е, о и

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

1. Ортогонализация векторного базиса 1.1. Точная ортогонализация

Обозначим Sp(a,b) пространство полиномов порядка не выше p, определенных на отрезке [a,b], со скалярным произведением

b

(u,v) = u(x) • v(x)dx. (2)

Пусть на отрезке [0,1] заданы точки Qpx = (жо, ...,жр} такие, что 0 = жо < ж1 < ... < жр-1 < жр = 1. Тогда набор интерполяционных полиномов Лагранжа из Sp(0,1) задается соотношением

p ( — )

LP(x.Q)= П ■ (3)

Определим параллелепипеидальный конечный элемент как множество

E = ((ж,y,z) | жо < ж < Жр,уо < У < Ур.^о < z < Zp}. (4)

Аналогичным образом введем Qy и Qp. Тогда простейший базис пространства функций из H(rot, E) на элементе (4) имеет вид

Nip)ijfc(ж,y,z) = Lp-1(x,Qx-1)Lp(y,Qp)LP(z,QP) ■ e*, N(p)ijfc(ж,y,z) = LP-1(y,QPp-1)Lp(z,QP)Lp(ж,Qx) ■ ey, (5)

NZp)ijk(ж,y,z) = Lp-1(z, QP-1)Lp(ж,Qx)Lk(y,Qp) ■ e^,

где г = 0,р — 1, = 0,р, еи — единичный орт, направленный вдоль оси координат и. Для упрощения дальнейшего изложения будем называть полиномы, определенные на оси координат, совпадающей с направлением базисной функции (т.е. ¿Р :(ж— в определении , ¿р-1(у,^у-1) — в определении и Ьр-1^,ОР-1) — в опреде-

лении ), нормальными, а все остальные тангенциальными.

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

Поскольку с ростом порядка базиса растет заполненность матрицы массы ненулевыми элементами, актуально построение такого базиса N1^, который, с одной стороны, обладал бы хорошими интерполяционными свойствами, с другой — минимизировал бы количество ненулевых элементов в матрице массы, увеличивал бы диагональное преобладание и снижал бы ее обусловленность. Поэтому ниже рассмотрим задачу выбора полиномов Ьр(ж, более подробно.

Можно заметить, что для конструкции базиса не имеет значения конкретный вид ¿к. Это позволяет варьировать их коэффициенты с целью получения базисных функций, обладающих свойствами, лучшими с точки зрения конечно-элементной дискретизации. В частности, для получения матриц массы с малым числом обусловленности используются полиномы, которые строятся неравномерным распределением узлов а = ж о < ж1 < ... < жр-1 < жр = Ь. Из теории аппроксимации известно [10, 11], что минимизация константы Лебега

р

А(д)= шах V ¿Р(ж ,д) (6)

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

= ес8[(2г + 1)п/(2р + 2)] 1 = _1 ж = — ес8[п/(2р + 2)] +2,г = 0,РГ (7)

Выбор полиномов Лежандра в качестве нормальных полиномов оптимален с точки зрения структуры матрицы массы. Полиномы Лежандра, помимо ортогональности, позволяют построить иерархический базис. На отрезке [0,1] они принимают вид

Тр(х) = ^/2^+1 ¡р-1(2х - 1), (8)

где ¡р(х) определяется на отрезке [-1,1] стандартным образом [12]:

10(х) = 1, ¡1(х) = х,

¡р+1(х) = 2р+1 х1р(х) -^—¡р-1(х). р + 1 р + 1

(9)

В этом случае векторный базис вводится как

(х,у,х) = ¡г-1(х)ЬРр(у,Я)Ьрк• ех,

(х,у,г) = ¡*-1(уЩ(г,Я)Ькк(х,Я) • еу, (10)

(х,у,г) = ¡г-1(г)Ьр(х,Я)ЬРк(у, Я) • ег.

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

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

Окончательно, необходимо построить две системы ортогональных полиномов на [0, 1]:

1) система нормальных ортогональных полиномов без ограничений (в качестве примера уже приведены полиномы Лежандра);

2) система тангенциальных ортогональных полиномов, обеспечивающая непрерывность интерполируемой функции в точках 0 и 1 .

Во втором случае требуется найти полиномы Ьг £ Яр(0,1), г = 0,р, такие, что

Ьо(0) = 1,Ьо(1) = 0,

Ьг(0) = Ьг(1) = 0,г = 1,р - 1, (11)

Ьр(0) = 0,Ьр(1) = 1.

Для этого Ьг, г = 0,р, следует искать в виде

Ь0(х) = (1 - х)М0(х),

Ьг(х) = х(1 - х)Мг(х), г =1,р - 1, (12)

Ьр(х) = хМр(х).

Здесь М0,Мр £ Яр-1(0,1), Мг £ Яр-2(0,1), г = 1,р - 1, подлежат определению. Имеет место следующее утверждение.

Утверждение 1. Ортогональная, т.е. обладающая свойством Уг,] = 0,р,

(Ьг,Ьз ) = 0, г = 3, (13)

(Ьг,Ь3 ) = 0, г = з, (13)

система полиномов вида (12) не существует. При этом существует система (12), обладающая свойством

(Ьг,Ьг) = 0, г = 0,р, _

I г = 1 р - 1

(14)

(Ьо, Ьг) = (Ьр, Ьг) = 0, г =1, р - 1,

(Ьг,Ьз) = 0 г = 3, = 1,р - 1

(Ьо,Ьр) = 0.

Иными словами, можно построить систему (12), в которой ортогональны все полиномы, кроме Ьо и Ьр (или любой другой пары).

Доказательство. Невозможность существования доказывается от противного. Пусть построена ортогональная система Ь0,..., Ьр-1. Будем искать Ьр в виде

/р-1

Ьр(х) = х I ^^ а3х3 \]=0

Тогда для определения Ьр получаем р уравнений вида (Ьг, Ьр) = 0 или У г = 0,р - 1:

р-1

(Ьг,х ) = 0. (15)

з=о

Матрица М с коэффициентами М^ = (Ьг,х3) — квадратная размерности р. Так как Ьг линейно независимы, х3 линейно независимы, строки М также линейно независимы. Следовательно, СЛАУ (15) с нулевой правой частью имеет только тривиальное решение.

Конструирование системы, удовлетворяющей условиям (14), состоит из следующих шагов:

1) строится ортогональная система Ь1,..Ьр-1 с помощью стандартной процедуры Грама—Шмидта. Поскольку линейная комбинация финитных функций финитна, то достаточно в качестве начальной системы взять любые линейно независимые полиномы вида х(1 - х)М3(х), М3 £ Бр-2(0,1), ] = 1,р - 1;

2) строится Ь0:

Шо = (1 - х)Мо(х), Мо £ Яр-1(0,1), Р-1

Яо = Шо -£ (Ьг,Шо)Шо, (16)

г=1

Ьо = Яо/Яо(0);

3) строится Ьр:

Шр = хМр(х), Мр £ Яр-1 (0,1), Р-1

Яр = Шр -Е (Ьг,Шр)Шр, (17)

г=1

Ьр = Яр/Яр(1). 1.2. Приближенная ортогонализация

Полиномы вг £ Бр(а,Ь), г = 1,к, называются ортогональными на системе точек г, £ [а,Ь], д = 1,т, если Уг = ], г,] = 1, к, д = 1,т:

вг(гя) • в3(г,) = 0. (18)

Идея построения приближенно ортогонального базиса состоит в замене скалярного произведения (2) формой

N

и(жг) ■ 1>(жг). (19)

г=0

Здесь точки жг и веса берутся из квадратур Гаусса или Гаусса—Лобатто [7]. Кроме того, полиномы ¿г(ж) £ Яр(а,6) из определения (5) строятся ортогональными на системе точек жг, в этом случае < >р+1 = 0, г = . Если порядок полиномов доста-

точно высок, то форма (19) достаточно точно аппроксимирует скалярное произведение (2). Для тангенциальных полиномов система точек жг должна включать точки 0 и 1, поэтому необходимо использовать квадратуры Гаусса—Лобатто [13]. Использование точек квадратуры Гаусса для построения нормальных полиномов целесообразно только для случая неортогональных сеток, в случае же ортогональной сетки можно использовать полиномы ортогональные в скалярном произведении (2).

2. Вычислительный эксперимент

Цель эксперимента — сравнение интерполяционных свойств различных базисов, а также вычислительных затрат при их использовании. Введем обозначения: Ир — стандартный базис, в котором в качестве нормальных и тангенциальных полиномов используются полиномы Лагранжа с корнями в узлах, определенных соотношением (7); ИРГ — базис, в котором полиномы Лежандра применяются в качестве нормальных полиномов и полиномы, определенные уравнениями (14), — в качестве тангенциальных; Ирг — базис, в котором нормальные полиномы — полиномы Лежандра, а тангенциальные полиномы — это полиномы, ортогональные на системе точек квадратуры Гаусса—Лобатто.

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

8т(жуг + вт(ж) есв(г)) есв(жуг — зт(у) есв(ж))

f (ж,У,г)

(20)

на единственном конечном элементе [0,1] х [0,1] х [0,1]. В такой задаче требуется решить СЛАУ Мж = г, где = (£, N1). Для данной симметричной положительно определенной матрицы А через £ (А) = тах( ) обозначим ее число обусловленности (здесь

Атт (А)

Атт(А) и Атах(А) — наименьшее и наибольшее собственные числа), через пг(А) — число ненулевых внедиагональных элементов матрицы А. Введем также обозначения: Мр — матрица массы для базиса Ир, Мрг — матрица массы для базиса ИРГ, М^ — диагональная матрица массы для базиса Ирг, ер, ерг и ерг — соответствующие базисам ошибки аппроксимации. В табл. 1 приведены значения числа ненулевых элементов в матрицах Мр и Мрг, числа обусловленности и нормы ошибок для базисов порядков р = 1, 6.

В табл. 2 приведены те же величины для разбиения области [0,1] х [0,1] х [0,1] на восемь равных кубических конечных элементов. Разбиение формируется делением каждой из сторон пополам. Для исследуемых базисов имеет место априорная оценка погрешности [3, 4]

||и — При||н (rot.fi) < СйР||и||(ЯР+1(П))3,

Таблица 1. Сравнение базисов для одного элемента

р N их(М Р) их(Мрг) £(М Р) £(М0Р) 1|еРЬ II Р 11 11 еог11Ь2 II V II ||Ь2

1 12 36 36 8.99 8.99 0.267318 0.267318 —

2 54 918 96 137.7 143.4 0.0326764 0.0326764 0.0332324

3 144 6768 180 261 394.7 0.00553697 0.00553697 0.00553939

4 300 29700 288 394.3 867 0.000851805 0.000851805 0.000853101

5 540 96660 420 703.5 1440 9.9662е-05 9.9662е-05 9.96709е-05

6 882 258426 576 1148 2401 1.27609е-05 1.27609е-05 1.27661е-05

Та блица 2. Сравнение базисов для восьми элементов

Р N их(М р) иг(Мрг) £(м р) е(мр) 1К|к2 II р II 1 |еог1 |Ь2 и р и ||ей1 ||Ь2

1 54 240 240 13.8 13.8 0.133513 0.133513 —

2 300 6636 672 65.4 91.4 0.00839971 0.00839971 0.00853654

3 882 51012 1296 182 240 0.000731928 0.000731928 0.000732591

4 1944 228552 2112 379.5 600 5.34082е-05 5.34082е-05 5.34827е-05

5 3630 752520 3120 684 1260 3.2368е-06 3.2368е-06 3.23774е-06

6 6084 2026260 4320 1125 2352 2.10231е-07 2.10231е-07 2.10297е-07

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

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

дН д2Н го1 го1 И + +

д* д*2

Н х п|дп = 0,

Г01/п

(21)

в области О = [-50, 50] х [-50, 50] х [-50, 50] (рис. 1). Источник поля — тонкая прямоугольная петля на границе областей О1 и О2. При этом = £(О2) = £0, = 0, а(О2) = 1. Задача в такой постановке описывает довольно широкий класс проблем, возникающих в геоэлектрике.

Сила тока задается функцией I(*), схематично представленной на рис. 2:

I (*)

/(*),* е [0,т], 1о,* е (Т, 4Т),

е [4Т,5Т], 0,* > 5Т,

(22)

где /(*) — монотонно возрастающая функция, /(0) = 0, /(Т) = 10, $(*) убывающая функция, $(4Т) = 10, ^(5Т) = 0.

После конечно-элементной дискретизации приходим к уравнению

«Н + м • ^Н + ме £ = г,

монотонно

где

я.

/ го1^ ■ гotNj ^О,

М?

г?

М? = / ■ N. ¿О,

п

= / ■ N. ^О,

п

/ гotJ0 ■ N. ^О.

п

/ /

Источник

п}/ // / /

п2 /

Рис. 1. Область моделирования О

Рис. 2. Зависимость силы тока от времени

С учетом того, что

у гоЛо • Ыг ¿0 = у Ло • го1Ыг ¿0 ^ (Ло х Ыг) ¿0

п п п

и

J (Ло х Ыг) ¿0 = ! (Ло х Ыг) • п ¿Б = 0,

П дП

вычисление {г для случая бесконечно тонкого проводника сводится к простому одномерному интегралу.

Явная схема по времени для уравнения (23) имеет вид

ЯИ +

1

АЬ

М * (Ига+1 - Ип) +

п+1

+М е

2

И

п+1

2

Ип +

^Аьп+1$п Аьп+1Аьп

где АЬП = Ьп - Ьп-1, в,п = АЬП + АЬп+1, или

И

п- 1

(24)

х

И

^ + -Я +

п+1 _

АЬ

1 М * + -Ме

п+1

АЬп+181

-1

х

АЬ

-М * +

2

п+1

АЬп+1Аьп

Ме Ип -

еттп-1

МеИ

(25)

Для решения задачи (23) будем использовать схему (25) в сочетании с базисом с неполной ортогонализацией и базисом, ортогональным на системе точек Гаусса—Ло-батто. Целью сравнения является установить:

— затраты памяти (число ненулевых элементов в матрицах массы и жесткости);

— возможность точного обращения матрицы массы и временные затраты на него. При дискретизации (23) используем сетку 10 х 10 х 10 элементов. Будем обозначать

через ввк(Л) — максимальное число элементов, которое может содержаться в матрице А-1. Фактически, это число равно числу элементов Л при переводе ее в профильный формат (более распространенное название формата — ББК [14]). Расчеты для стандартного базиса не проводились, поскольку его использование заведомо более затратно

п

1

2

Таблица 3. Сравнение характеристик матриц для различных базисов

p N nz(M'r) nz(R'r) nz(RPQl) ssk(M'r) t or ,c tgl ,c

2 21660 51840 861 600 476 064 15 898 060 2151 51

3 75 690 106 560 4360 680 2825 064 12 6336 627 8753 169

4 182 520 180 480 14 563 200 10 348 416 544 551216 35 352 461

5 360150 273 600 38 220 600 28 841400 1684 461475 114 796 1091

6 626 580 385 920 85 554 720 67339 296 4229 360172 279 238 2324

и не дает никакого выигрыша в точности. В табл. 3 приведены данные о размерностях и числе ненулевых элементов в матрицах массы и жесткости, а также время решения для 7500 итераций по времени с шагом 10-10 с. Значение T полагалось равным 10-7, a = 1. При использовании базиса B'r для решения СЛАУ с матрицей M' использовался метод сопряженных градиентов с требованием уменьшения невязки в 1016 раз и ограничением на число итераций, равным 2000. При реализации использовалось распараллеливание основных операций (умножения матрицы на вектор, линейных комбинаций векторов, скалярного произведения и пр.) с помощью технологии OpenMP. Результаты измерений времени приведены для тестов на машине с четырехъядерным процессором Intel Core 2 Quad 6600 с распараллеливанием на четыре потока.

Из результатов вычислений можно сделать следующие выводы:

— точное обращение матрицы массы нецелесообразно, поскольку число элементов в матрице ssk(Mo') оказывается чрезвычайно высоким несмотря на малое исходное число элементов в матрице M (действительно, полное LLT- или LDLT-разложение матрицы, содержащей примерно 0.5 • 109 элементов, потребует 1013 — 1014 сложений и умножений, а обратная подстановка — 0.5 • 109 сложений и умножений, что существенно сложнее, чем даже многократное решение СЛАУ с матрицей M);

— решение задачи с использованием базиса B'r приблизительно в 50-120 раз более затратно по времени вычислений и примерно на 30-50 % — по объему памяти. Естественно, что подобные затраты связаны с требованием достаточно качественного решения СЛАУ на каждом шаге. В большинстве реальных случаев достаточно решения СЛАУ с требованием уменьшения нормы невязки в 108 — 1010 раз, кроме того, на скорость решения СЛАУ большое влияние оказывает выбор начального приближения (которое можно находить, например, экстраполяцией), а также выбор предобусловлива-теля. Исследование этого вопроса, предположительно, позволит значительно сократить время решения СЛАУ с матрицей M' и сделать его сравнимым с временем умножения матрицы R на вектор.

На рис. 3-7 показаны графики зависимостей Hz (t) в точке (5,0, 0) на дневной поверхности. Как видно, решения, полученные на базисах Bgr и B^, значительно различаются. Кроме того, качественная картина поведения поля для этих двух базисов не соответствует физике процесса: на отрезке t £ [0,10-6] величина |Hz(t)| должна монотонно возрастать, чего не видно на рис. 3. Из этого можно сделать вывод, что использовать базис B2gl нецелесообразно в силу низкой точности. Очевидно также, что применение базисов пятого и шестого порядков не дает особых преимуществ, однако чрезвычайно усложняет вычислительную процедуру. Можно, таким образом, заключить, что по совокупности точности и вычислительных затрат для задач подобного класса лучше всего использовать базис B4i (из рис. 4 видно, что в случае Bgi ошибка, обусловленная конструкцией базиса, все еще достаточно велика) и при необходимости измельчать

Рис. 3. Н(£), £ = 107£, в точке М для базисов Рис. 4. Н(£), £ = 107£, в точке М для базисов

Бог (график 51) и Б^ (график 52)

во

(график 51) и Б^ (график 52)

Рис. 5. Н(£), £ = 107£, в точке М для базисов Рис. 6. Н(£), £ = 107^, в точке М для базисов

Бог (график 51) и Б^ (график 52)

во

(график ^1) и Б^ (график 52)

/ / / / / / / !

\ \ \ \ \ \

! I \ I \ \

Рис. 7. Н(£), £ = 107£ в точке М для базисов Рис. 8. Н в сечении плоскостью у = 0 для

Бог (график 51) и Б^ (график 52)

базиса Б?г

пространственную сетку. На рис. 8 приведен вид решения H в сечении плоскостью y = 0 для базиса B'r и p = 3, полученный аналогичным методом на сетке 30 х 30 х 30 элементов.

Заключение

В работе введен и исследован базис пространства H(rot, П), который строится с использованием частично ортогональных полиномов, что позволяет значительно уменьшить вычислительные затраты на решение уравнений Максвелла на ортогональных сетках. Проведено сравнение с базисом, предложенным в [7] для волнового уравнения, который был вполне успешно применен к задаче в области, содержащей a = 0.

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

[1] BossAviT A. Computational electromagnetism. Variational formulation, complementary, edge elements. San Diego: Acad. Press, 1998.

[2] Hiptmair R. Finite elements in computational electromagnetism // Acta Numera. 2002. N 11. P. 237-339.

[3] Nedelec J.C. Mixed finite elements in R3 // Numer. Meth. 1980. Vol. 35. P. 315-341.

[4] Nedelec J.C. A new family of mixed finite elements in R3 // Numer. Meth. 1986. Vol. 50. P. 57-81.

[5] Graglia R.D., Wilton D.R., Peterson A.F. Higher order interpolatory vector bases for computational electromagnetics // IEEE Trans. Antennas and Propagation. 1997. Vol 45, N 3. P. 329-342.

[6] Rodrigue G., White D. A vector finite element time-domain method for solving Maxwell's equations on unstructured hexahedral grids // SIAM J. Sci. Comp. 2001. Vol. 23, N 3. P. 683-706.

[7] Fisher A., Rieben R.N., Rodrigue G.H., White D.A. A generalized mass lumping technique for vector finite-element solution of the time-dependent Maxwell equations // IEEE Trans. on Antennas and Propagation. 2005. Vol. 53, N 9. P. 2900-2910.

[8] Jensen M.S. High convergence order finite elements with lumped mass matrix // Int. J. Numer. Meth. Eng. 1996. Vol. 39, N 2. P. 1879-1888.

[9] He B., Teixeira F.L. Sparse and explicit FETD via approximate inverse hodge (mass) matrix // IEEE Microwave and Wireless Components Letters. 2006. Vol. 16, N 6. P. 348-350.

[10] Brutman L. Lebesgue functions for polynomial interpolation — a survey // Annals of Numer. Math. 1997. Vol. 4. P. 694-704.

[11] Rieben R., White D., Rodrigue G. Generalized high order interpolatory 1-form bases for computational electromagnetics // Proceedings of the 2002 IEEE International Antennas and Propagation Symposium. San Antonio, Texas. 2002. Vol. 4. P. 686-689.

[12] Vosse F.N. van de, Minev P. Spectral element methods: theory and applications // Eindhoven University of Technology Report 96-W-001 ISBN 90-236-0318-5. 1996.

[13] Бахвалов Н.С. Численные методы. М.: Наука, 1973. Т. 1.

[14] Saad Y. Iterative Methods for Sparse Linear Systems. Boston: PWS Publishing Company, 1996.

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