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

Математическое моделирование

УДК 517.958, 530.145.6, 519.632.4 БГО: 10.22363/2312-9735-2017-25-3-217-233

Метод конечных элементов высокого порядка точности решения краевых задач для эллиптического уравнения в

частных производных А. А. Гусев

Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московская область, Россия, 141980

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

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

1. Введение

В работах [1,2] разработаны символьно-численные алгоритмы и программы решения краевых задач для системы обыкновенных дифференциальных уравнений второго порядка методом конечных элементов (МКЭ) высокого порядка точности с применением интерполяционных полиномов Эрмита (ИПЭ). Такая реализация МКЭ в отличие от традиционной с использованием интерполяционных полиномов Лагранжа (ИПЛ) обеспечивает непрерывность производных до заданного порядка приближенного решения не только на конечных элементах сетки, но и на границах конечных элементов [3], т.е. обеспечивает законы сохранения тока носителей заряда в квантово-размерных полупроводниковых системах или тока вероятности в кванто-вомеханических задачах рассеяния [4]. Вычислительная схема МКЭ решения с высокой точностью краевых задач на собственные значения для эллиптического уравнения в частных производных в ограниченной области двумерного пространства с применением ИПЛ дана в работе [5]. Обобщение вычислительной схемы МКЭ для решения эллиптических краевых задач с применением ИПЭ составляет цель данной работы.

Статья поступила в редакцию 31 марта 2017 г.

Автор благодарит О. Чулуунбаатара, С.И. Виницкого, В.Л. Дербова, В.П. Гердта, А. Гужджа и

Л.А. Севастьянова за сотрудничество.

Работа поддержана РФФИ (грант 16-01-00080).

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

Структура работы следующая. В разделе 2 дана постановка задачи. В разделе 3 дана формулировка вычислительной схемы МКЭ с применением ИПЭ нескольких переменных. В разделе 4 представлены результаты расчётов эталонной краевой задачи, демонстрирующие эффективность вычислительной схемы и порядок точности по числу конечных элементов разбиения области и числу кусочно-полиномиальных функций при различных ИПЭ двух переменных. В заключительном разделе даны перспективы развития и применения предложенных вычислительных схем.

2. Постановка задачи

Рассмотрим самосопряжённую краевую задачу для эллиптического дифференциального уравнения второго порядка:

(D - Е) Ф(х) = | -£ А9ij(z)+ V(z) - Е | Ф(х) = 0.

9o(z) i/-=1 dzi

(1)

Для коэффициентов главной части (1) выполняется условие равномерной эллиптичности в ограниченной области г=(г1,..., евклидова пространства т.е. существуют константы ^ > 0, и > 0 такие, что

d £

i,3 = l

i=l

Также предполагается, что g0(z) > 0, gji(z) = gij(z) и V(z) — функции, непрерывные вместе со своими обобщёнными производными до заданного порядка в области, z G Q = Q U дО с кусочно-полиномиальной границей S = дО, обеспечивают существование нетривиальных решений, подчинённых граничным условиям первого (I), второго (II) или третьего (III) рода [6,7]:

(I): Ф(г)|5 = 0,

(II) : 8*(Z)

(III) :

дпв д Ф(г

0,

дп

D

+а(в)Ф(г)|5 = 0,

дФ(г) дп

D

I лдф(г) (n,ei)gij(z)-

i, 3 = 1

(2)

dz.

где а(з) = а(з(х)) — известная функция на границе области з(г) € Б, дФт(¿)/дпв — производная по направлению конормали, п — внешняя нормаль к границе области

а

Б = дО, ег — единичный вектор вектора г = ^ пггг, ("п, п) — скалярное произведе-

г= 1

ние в

Для задачи дискретного спектра собственные функции Фт(г) из пространства

Соболева Н^^П), Фт(г) € Нсоответствующие собственным значениям энергии Е: Е1 ^ Е2 ^ ... ^ Ет ^ ... удовлетворяют условиям нормировки и ортогональности

{§т^)1<Фт'(г)) = !до(г)Фт(г)Фт'(г)йг = 5тт', йх = dzl (3)

п

Решение МКЭ краевых задач (1)—(3) сводится к нахождению стационарных точек вариационного функционала [7,8]

^(Фт, Ет, £) = I до(г)Фт(г) (Б — Ет) Ф(г)йг = п

= П(Фт,Ет, г) Фт(г)^¡^йз, (4)

где П(Ф

то Еm, г) — симметричный квадратичный функционал вида П(Фт,Ет, £) =

<

Поверхностный интеграл в (4) вычисляется с учётом (2).

Е 9ц( ¿) 9о(г)Фт(г)Фт(г)(У(г) — Ет)

г,] = 1 г °

йх.

(5)

3. Вычислительная схема МКЭ

В МКЭ область П = Пк(г) = и<^=1 Ач, заданная в виде многогранника, покрывается конечными элементами, в данном случае симплексами Ая с й + 1 вершинами ' = (¡г1, ,..., при г = 0,..., й. Каждое ребро симплекса Ад разбиваем на р равных частей и проводим семейства параллельных гиперплоскостей Н(I,к), нумеруя каждую целым числом от к = 0,..., р, начиная от соответствующей грани, например, как показано при й = 2 на рис. 1 (см. также [9, С. 220]). Уравнение гиперплоскости Н(I, к) : Н(г; £) — к/р =0, где Н(г; х) - линейная функция от г.

Узловые точки пересечения гиперплоскостей Аг нумеруем наборами целых чисел [п0,... ,пл], щ > 0, п0 + ... + па = р, где т, г = 0,1,... — номера гиперплоскостей, параллельных грани симплекса, не содержащей г-вершину ¡¿1 = (2ц,... '^о). Координаты ¡1Г = (£Г1,..., узловой точки Аг € Ая вычисляются по формуле

(&1,..., €гсд = (¡01,..., ¡оа)по/р + (¿11,..., ¿1а)п1/р + ... + (¿<11,..., '¡аа )па/р (6)

через координаты вершин ¿^ = ('1,..., ¿^а). Тогда ИПЛ рг(¿) равные единице в точке Аг с координатами = (£г1,..., характеризуемой числами [п0,п1,... ,по], и нулю в остальных точках £г', т.е. рг() = 5ГГ' имеют вид

<М.Ь |П П НЩ| г (7)

(Ь)

Рис. 1. Нумерация узлов Аг, г = 1,..., (р + 1)(р + 2)/2, наборами чисел [п0,п1,п2\ стандартного треугольного элемента Д для схемы с ИПЛ пятого порядка р' = р = 5 при d = 2; линии (пять пересекающихся прямых) — нули ИПЛ ^14(г') из (13), равного единице в точке, пронумерованной тройкой чисел [п0,п1,п2\ = [2, 2,1\; (Ь) изолинии ИПЛ р14(г')

Отметим, что построение ИПЭ ^'^(х) с фиксированными значениями функций {^г(&')} и производных и=£г'} в узлах уже при d = 2 приводит к гро-

моздким выражениям, затратным для использования в МКЭ на неравномерных сетках.

Экономичная реализация, принятая в МКЭ, следующая:

— расчёты проводятся в локальных координатах х' , в которых координаты вершин симплекса А следующие: х'^ = ..., ^ = &зк,

— ИПЭ в исходных координатах х ищется в виде линейных комбинаций полиномов в локальных координатах х' , при этом переход в исходные координаты выполняется только на этапе численного решения конкретной краевой задачи

(1)-(4),

— вычисление интегралов МКЭ выполняется в локальных координатах.

Построим ИПЭ на произвольном ¿-мерном симплексе Ад с <1 +1 вершинами

= (2ц..., ¿го), г = 0,... Для этого введём локальную систему координат х' = (х\, х'2,..., х'л) € в которой координаты вершин симплекса следующие: х\ = (х[к = 5ik ,к = 1,... ,д). Связь между координатами даётся формулой:

= ¿0г + - ¿0г , г = 1,...,^

3 = 1

(8)

или в матричном виде

г02

+3

4

, з=

\zodJ \хУ

(Хц-Х01 г21-Х01 Х!2-Х02 ¿22-¿02

г<ц-г01\ ¿й2—%02

(9)

г

Обратное преобразование имеет вид

1 1 /¿оЛ

^-1 —J-1 ¿02

\Z0dJ

(10)

Отсюда следует связь между операторами дифференцирования

/_дД д*1 \ д дг2 =.]т /_а_\ дх1 _а_ дг2 /_а_\ дх1 _а_ дг2 =0-1)Т /0\ д*1 д д4

\ \ш2/ \~dki)

(11)

Формула (11) применяется для вычисления ИПЭ х')=(ф"(г')^3(г')] из (21), удовлетворяющих условиям (14), (18), (19) из следующего раздела, с фиксированными производными до заданного порядка в узлах £г'. При этом производные по нормали к границе элемента в исходной системе координат в общем случае не являются таковыми в локальных координатах х'. При построении ИПЭ в локальных координатах 4 требуется пересчёт фиксированных производных в узлах £г' элемента Ая в узлы £'г' элемента А, используя матрицы J-1, заданные громоздкими выражениями. Поэтому требуемый пересчёт производных выполняется, используя соотношения (8)—(11), для каждого конечного элемента Ая на этапе формирования

базиса ИПЭ (у" (на конечном элементе А, который выполняется численно по аналитическим формулам, представленным в следующем разделе.

Интегралы, входящие в вариационный функционал (4) на области ^ь(^) = Ад, выражаем через интегралы, вычисленные на элементе Ад, пересчитанные в локальные координаты 4 на элементе А,

I до(г)<р?(¿)укг, (г)и= ^ д0(г(г'))укг(г')укг, (г')и(г(г'))йг',

Д0 А

9*1*2 (Х)~

д^(г) дуУ (г) дг81 д хВ2

йг =

(12)

^ Е ^ 1) ¿1*1 ^ 1) ^2]

tl, Ь2 = 1 а

д 4

д

2

где J=<iet J>0 определитель матрицы J, у"(г')=(фг(г')^3(г')] из (21), йг' =

йх 1 ... й .

3.1. Интерполяционные полиномы Лагранжа

В локальных координатах ИПЛ уг(4), равные единице в узловой точке , характеризуемой числами [п0,п1,..., по\, и нулю в остальных узловых точках £'г', т.е., У г (Сг') = $гг' определяются по формуле (7) при Н (0; г')=1 —г[ — ... — Н (г, х')= г=1,...

у м= (п п 4-%1 (п1-г1-у I ■ (13)

\1=1^0щ/р-пзр) \п0=0 п0/р-щ/р

Приравнивание числителей в (13) нулю даёт семейства уравнений прямых, направленных «горизонтально», «вертикально» и «наклонно» в локальной системе координат элемента А, которая связана аффинным преобразованием (8) с семействами «наклонно» расположенных прямых элемента Ая. На рис. 1 дан пример, иллюстрирующий построение ИПЛ при й =2, г, г' = 1,..., (р +1)(р +2)/2, р = 5 на элементе А в виде прямоугольного треугольника с вершинами х'0 = (х'01, х'02) = (0, 0), #1 = (2'и, = (1,0), ¿2 = (¿21, ¿22) = (0,1).

Кусочно-полиномиальные функции Рц(г), формирующие конечноэлементный базис {Р]т(х)}^=1, которые строятся путём сшивки ИПЛ (х) из (7), получаемых из (13) с помощью преобразования (10), на конечных элементах Ад:

Рг(х) = {^¡(х), А1 € Ад; 0, А1 € Ад} ,

являются непрерывными, но их производные терпят разрывы на границах элементов Аа.

3.2. Интерполяционные полиномы Эрмита

Построим ИПЭ порядка , сшивкой которых можно получить кусочно-полиномиальные функции (25), сохраняющие непрерывность производных до заданного порядка К.

3.2.1. Вспомогательные полиномы (ВП1)

Для построения ИПЭ в локальных координатах введём набор вспомогательных полиномов (ВП1)

,,К\...К<1(С' \_

(Ргг'"Ка () = 5гг' 5К10 . . . ,

дх 111 ...д х'/*

=5ГГ'5К1^1 ... ,

*>=К' (14)

0 < К1 + К2 + ... + ка < ктах-1, 0 < + Ц2 + ... + Ца < ктах-1.

Здесь в узловых точках , определённых согласно (6), в отличие от ИПЛ заданы значения не только функций, но и их производных до порядка ктах-1. ВП1 даются выражениями

у) = Юг у) £ (,'1 - ^ х ... х (^ - ^у. ,

п, (,)= ( А "л ^ г-</р)^ \ (Ъ1 (1-г 1- ... - 4-п0/рГ- I (15) {11п\=0 (п*/р-</р)к™) Ц=0 (п0/р-п'0/рг-

(С) =1

где коэффициенты вычисляются из рекуррентных соотношений, полу-

ченных в результате подстановки (15) в условия (14),

' 0, ¡л1+ ■■■ < ... +ка, (ц1= (к1,...,ка); а г

Ц—!, (¡1,..., ¡а) = (к1,...,ка);

£ (й ^ - "<)')

Уг г ,

(16)

¡1+ ... > К1+ ... +ка;

2-к* (г') =

1 дК1К2...Кл,Шг (^ шг(г') дх[к1д.. .дг'а

При д > 1 и ктах > 1, число

ИПЭ степени р' и кратности узлов кт

меньше числа полиномов N1р , формирующих базис в пространстве полиномов степени р' (например, ИПЛ из (13)), т.е. полиномы, удовлетворяющие (14), определены неоднозначно.

3.2.2. Вспомогательные полиномы (ВП2 и ВП3)

Для однозначного определения полиномиального базиса введём К = N1р — вспомогательных полиномов Qs(z) двух типов: ВП2 и ВП3, линейно-независимых от ВП1 из (15) и удовлетворяющих условиям в узловых точках £Г' ВП1:

Qs^Г') 0 д¿1гдг'12 ...дх'

1в

К,

(17)

0 < К1 + К2 + ... + Ка < Ктах — 1, 0 < ¡1 + ¡2 + ... + ¡а < Кт&х — 1.

Отметим, что для обеспечения непрерывности производных часть полиномов, называемых ВП2, должны удовлетворять условию

дkQs(z')

дп\з)

=5ss', в, в'=1,... ,Т1(к ), к = к(в'),

(18)

г' =р

где г]'^ = (т]'^ 1,..., т/'^а) — выбранные точки, лежащие на гранях всевозможных размерностей (от 1 до д — 1) д-мерного симплекса А и не совпадающие с узловыми точками ИПЭ , в которых выполняется (14), д/дщ(8) — производная по направлению вектора щ, направленного по нормали к соответствующей г-й грани д-мерного симплекса Ая в точке г]' в исходной системе координат, который пересчитывается в точку г/', грани симплекса А в локальной системе координат, используя соотношения (8)—(11), например, при д =2 см. (23). Вычисляя число Т1(к) независимых параметров, требуемых для обеспечения непрерывности производных до порядка к, определяем его максимальное значение к , которое можно получить для схем с заданными р и ктах, и, соответственно, дополнительных условий (18).

Остаётся Т'=К—Т1(к') независимых параметров и, соответственно, добавляется Т' дополнительных условий, необходимых для однозначного определения полиномов, называемых ВП3,

Qs (О )=б''', 8,8'=Т1(к' )+1,...,К,

(19)

г

1

где

Сз> (Сз >и ■ ■ ■, Св'а) е А выбранные точки, принадлежащие симплексу без гра,-ницы, но не совпадающие с узловыми точками ВП1 . Вспомогательные полиномы ВП2 даются выражением

(,&) = ^ ... х'ак* (1 - 4 - ■■■ - х'а)к° £ ... гТ, (20)

31,-; М

где к = 1, если точка г]8, в которой заданы дополнительные условия (18), лежит на соответствующей грани симплекса А, т.е. Н(1, г]8) = 0, и кг = К, если Н(1, г]8) = 0. Вспомогательные полиномы ВП3 даются выражением (20) при кг = К. Коэффициенты определяются из системы линейных уравнений, полученных в результате подстановки (20) в условия (17)—(19).

В результате получаем требуемый набор базисных ИПЭ

= (21) составленных из полиномов (8(х') типа ВП2 и ВП3 и полиномов ф^(х') типа ВП1:

дк^г(г')

К I ^ ^г

г;(х') = <£г;(х') - ^ СК,Г;вЯв(^), СК;г;8 = ^ дпк{8)

(22)

I ~ ' —''5

8=1

я8(х') е ВП2,

г' =Ч'В

Сз), Яэ(^) е ВП3.

3.2.3. ИПЭ при й = 2

При й = 2 степень р' полинома по тангенциальной переменной £ на границе треугольника Ая совпадает со степенью полинома двух переменных, и для его однозначного определения требуется р' + 1 параметр. Производная порядка К по нормальной переменной к границе будет полином степени р' -К, и для её однозначного определения потребуется р'-К + 1 параметров. Однако, она определена только р'-К(р +1) параметрами: смешанными производными фиксированного порядка К по нормальной переменной и порядка от 0 до Ктах -К - 1 по тангенциальной переменной, т.е.

дК> дК>+1 д^тах

дпК>' дпК>дЬ1 ' дпК>дЬКтах-К>-1

Таким образом, с учётом того, что у треугольника три стороны, для однозначного определения производных до порядка К на границе необходимо Т1 (К) = 3р + ... + Зп'р = 3рп'(К + 1)/2 параметров и, соответственно, дополнительных условий (18).

Например, для р = 1 и Ктах = 4 имеется К = 6 дополнительных условий для определения ВП2 и ВП3. Степень р' = 7 полинома по тангенциальной переменной £ на границе треугольника совпадает со степенью полинома двух переменных, и для его однозначного определения требуется р' + 1 = 8 параметров. Производная первого порядка К = 1 по нормальной переменной к границе будет полином степени р' -К = 6, и для её однозначного определения потребуется р' -К + 1 = 7 параметров. Однако она определена только р' - К(р+ 1) = 6 параметрами: смешанными производными д/дп, д2/(дпд1), д3/(дпд12), заданными в двух вершинах. Недостающий параметр можно определить, задавая производную по направлению, непараллельном границе треугольника в одной из точек на его стороне (например, посередине). Таким образом, для р = 1 и Ктах = 4 можно построить ИПЭ с фиксированными значениями первой производной на границе треугольника, при этом остаётся 6 - 3 = 3 свободных параметров.

Производная второго порядка К = 2 по нормальной переменной к границе — полином степени р' - К = 5, и для её однозначного определения потребуется

р' - К + 1 = 6 параметров. Однако она определена только р' - К(р + 1) = 4 параметрами: смешанными производными д2/дп2 и д3/(дп2дЬ), заданными в двух вершинах треугольника. Таким образом, для однозначного определения второй производной потребуется 6 параметров. Это означает, что данным алгоритмом для р = 1 и Ктах = 4 нельзя построить схему МКЭ с непрерывной второй производной. В этом случае следует использовать схему с Ктах > 4, например, обозначенную [152] из табл. 1 и рис. 2. Тогда оставшиеся три свободных параметра используются для построения ВП3. Заметим, что возможно построение схем, обеспечивающих непрерывность вторых производных на некоторых границах конечных элементов. Этот случай в данной работе не рассматривается.

При й =2 производные д/дп^ по направлению п^, перпендикулярном к соответствующей грани = 0, 1, 2 в исходной системе координат, даются через частные производные д/дх,, ] = 1, 2 в локальной системе координат А, используя (8)-(11), выражениями

д д д д д д дЩ+Ы' '=1'2' (23)

где Д,- = ¿0, ¡¿1, ¿2) - функции от координат вершин ¿0, ¿1, ¿2 треугольника А д в исходной системе координат

¡11 = 3 — 1 (( ¿21 ¿01)2 + (Х22-¿02)2)1/2,

Г = ( ¿12-¿02)( ¿22- ¿02) + ( ¿21 - ¿01)( ¿ц- ¿01)

Н2= 3((Х21-Х01)2 + (Х22-Х02)2)1/2 ' /22 = 3-1((Х11-¿01)2 + (Х12-¿02)2)1/2,

, = (¿12-¿02)( ¿22-¿02) + ( ¿21 - ¿01)( ¿11- ¿01) (24)

/21= 3((Х11-Х01)2+(Х12-Х02)2)1/2 1 ( )

¡02 =

/о1 = -(23 )-1((Х22 - ¿12)2 + (Х11-¿21)2)1/2, (¿11- ¿01 )2 + ( ¿12 - ¿02)2-( ¿22-¿02)2-( ¿21-Х01)2

23 ((¿22-¿2l)2 + (¿11-¿12)2)1/2 3 =(¿11- ¿01)^22-¿02)-(¿12-¿02)^21-¿01).

Схематическое изображение реализации условий (14), (17), (18) и (19), с помощью которых построены базисные ИПЭ, показано при й = 2 на рис. 2. Характеристики полиномиального базиса из ИПЭ на элементе А при й = 2 приведены в табл. 1. Выражения для ИПЭ (21) до р' = 9 порядка при й =2 представлены в Приложении.

3.2.4. Кусочно-полиномиальные функции

Кусочно-полиномиальные функции Рг(х) с непрерывными производными до порядка К строятся путём сшивки полиномов ^(£) = {фг(х),(8(х)} из (21), получаемых алгоритмом из Приложения, на конечных элементах Ад е г) = и^=1 Ая:

Р1>(г) = {(г),Акп е Ая;0,АЩ1) £ Ач) , (25)

где знак «-» может появиться только для ВП2, когда нужно сшить нормальные производные нечётного порядка.

Рис. 2. Схематическое изображение условий на элементе Дд для построения базиса ИПЭ [рктахк']: [131], [141], [231], [152]. Квадраты: точки , в которых фиксируются значения функций и их производных согласно условиям (14), (17); сплошные (прерывистые) стрелки: с началом в точках ^, в которых фиксируются значения первой (второй) производной по направлению нормали в исходных координатах, согласно условию (18), соответственно; кружки: точки , в которых фиксируются значения функций согласно условию (19)

Таблица 1

Характеристики базиса ИПЭ (21) при <1 = 2

[рктахк ] [131] [141] [231] [152]

р' Ктах(р + 1) — 1 5 7 8 9

^тахУ (р + 1)(р + 2)ктах (ктах + 1)/4 18 30 36 45

N1р/ (р' + 1)(р' + 2)/2 21 36 45 55

К р(р + 1)ктах (Ктах — 1) / 4 3 6 9 10

Т1(1) 3р 3 3 6 3

Щ2) 9р 9 9 18 9

N (ВП1) ^ „> ^тахУ 18 30 36 45

N (ВП2) Т1(К) 3 3 6 9

N (ВП3) К — Т1(К) 0 3 3 1

Ограничение на порядок производных К : 3рк' (К + 1)/2 ^ К.

Разложение искомого решения Фт(г) по базису кусочно-полиномиальных функций Ру(г), Ф^^) = ^=1Р1>(г)Ф$т и подстановка его в вариационный функционал (4) приводит к обобщённой алгебраической задаче на собственные значения, (А — ВЕ1 )§1п = 0 , которая решается стандартным методом (см., например, [8]). Элементы симметричных матриц жёсткости А и масс В содержат интегралы типа

(4), которые вычисляются на элементах в области Ая £ = и^=1 Ад, пересчи-

танные в локальные координаты на элементе А.

Теоретические оценки погрешностей между точным Ет, Фт(г) £ и при-

ближенным Е1, Ф^Кг) £ Н, решениями следующие [10,11]:

|ЕГО — Е1 < С1 Ь2р', \\Фт(г) — Фнт(г)\\0} < С2ЬР'+1, (26)

где

Ш*)\\0 = \ д0(г)йгЩг)Фг(г), п

Н — максимальный размер конечного элемента Ад, р' — порядок схемы МКЭ, т — номер собственного значения, С1 и С2 — коэффициенты, не зависящие от Н.

4. Результаты и обсуждение

В качестве примера рассмотрим решение задачи дискретного спектра (1)-(3) при й = 2, д0(г) = дц(х) = 1, и V(г) = 0 в области = и®=1Ач в виде рав-

ностороннего треугольника со стороной 4к/3 с граничными условиями второго типа (II), который разбит на ( = п2 равносторонних треугольников Ад со стороной Н = 4л/3п. Собственные значения этой задачи с вырожденным спектром — целые числа Ет = 0,1,1, 3,4,4, 7, 7,....

На рис. 3 представлены конечноэлементная сетка с ИПЛ восьмого порядка и профиль четвёртой собственной функции Ф\(х).

(b)

Рис. 3. Сетка на области = УДч треугольной мембраны,

составленная из треугольных элементов Дч; (Ь) профили четвёртой собственной функции Ф= 3, полученной на ИПЛ

порядка р' = р = 8

На рис. 4 показаны погрешности АЕт = Е1^ — Ет собственного значения E%(z) в зависимости от числа п (число элементов п2) и от длины N вектора Ф^ алгебраической задачи на собственные значения для схем МКЭ от пятого до девятого порядка точности, использующие ИПЛ, отмеченные метками [рктахк'] = [510],..., [910], и использующие ИПЭ, отмеченные метками [131], [141], [231] и [152] из табл. 1, сохраняющие непрерывность первой и второй производных приближённого решения соответственно.

Как видно из рис. 4, погрешности собственного значения АЕ1^) схем МКЭ одного порядка р' = Kmax(p+ 1) — 1 примерно одинаковые и соответствуют теоретическим оценкам (26), но в схемах МКЭ с ИПЭ, сохраняющих непрерывность первой и второй производных приближённого решения, используются матрицы меньшей размерности, соответствующие длине вектора N, в 1,5-2 раза меньшей, чем для схем с ИПЛ, сохраняющих на границах конечных элементов только непрерывность приближённого решения. Вычисления проводились на компьютере 2 x Xeon 3.2 GHz, 4

GB RAM, используя Intel Fortran 77 с четверной точностью real*16, с 32 значащими цифрами. Время счёта рассмотренных примеров не превышало 3 минут.

1Е-5

■ч*

ш

<1

1Е-10 1Е-15

1Е-20

—— [510]- — - [131] — [610] —■—[710]--«- - [141] ГЧ —— [810h-e-- [231] —— [910]-*- [152]

Ч\ ч

\ 4 \ ^f«

Ч Xvvcv

100 1000 п N

Рис. 4. Зависимости погрешности АЕ% собственного значения от числа элементов п и от длины вектора N

5. Заключение

Предложена вычислительная схема МКЭ высокого порядка точности решения задачи на собственные значения для эллиптического уравнения в частных производных в ограниченной области многомерного евклидова пространства, обеспечивающая непрерывность не только приближенного решения, но и его производных до заданного порядка. На примере точно-решаемой краевой задач для треугольной мембраны показано, что для достижения заданной точности приближённого решения, для схем с МКЭ с ИПЭ, обеспечивающих непрерывность первой и второй производных приближённого решения, используются матрицы меньшей размерности, соответствующие длине вектора N, в 1.5-2 раза меньшей, чем для схем с ИПЛ, сохраняющих на границах конечных элементов только непрерывность приближённого решения.

Отметим, что явное выражение для ИПЭ нескольких переменных, в отличие от ИПЛ, достаточно сложное, в локальных координатах 4 требуется пересчёт фиксированных производных в узлах хг' элемента Ая в узлы г'г' элемента А, используя

матрицы J-1. Поэтому требуемый пересчёт производных выполняется для каждого конечного элемента Ад на этапе формирования базиса ИПЭ {у^ (^')}^=1 на конечном элементе А, который выполняется численно по аналитическим формулам, аналогичным формулам при <1 =2 из раздела 3.4 и Приложения.

Вычислительные схемы МКЭ ориентированы на расчёты спектральных и оптических характеристик квантовых точек и других квантовомеханических систем. Реализация МКЭ с ИПЭ в пространстве <1 ^ 2 и других областях, отличных от многогранника, будет дана в последующих работах.

6. Приложение

Выражения для ИПЭ в исходных координатах, удовлетворяющие условиям (14), (18), (19), так, что для каждого полинома в правых частях этих условий единица появляется только один раз, громоздкие. Их вычисление удобно проводить по алгоритму:

1) по формулам (15) и (20) вычисляются ВП1 ^(г'), ВП2 и ВП3 ^(г');

2) используя формулу (22), вычисляем базис ИПЭ, состоящий из полиномов ф^(х') и Qs(z'), в локальной системе координат;

3) ВП1 <Рг(г) пересчитываются по формулам (11);

4) выполняется преобразование (10).

В табл. 2-5 приведены ВП1 ВП2 и ВП3 Ql|(г'), а соответствующие коэф-

фициенты ск;г;8 вычисляются по формулам (22). Обозначения: £г, г]3, Св — координаты узлов в которых правая часть (14), (18) или (19) равна единице, = 1 — — х2, определяются из формул (24), при этом аргументы у функций и штрихи у независимых переменных опущены.

Таблица 2

ИПЭ р =1, «тах = 3, к' = 1, р' = 5 (элемент Аргириса [3])

ВП1 : а = (0,1), & = (1, 0),6=(0, 0)

2,° 2 3 /о , °=х3(6х2-15х1+10) ^2 , 1=-Х3Х2(3Х1-4) ^2,°=-г13(^1-1)(3^1-4) °,2 3 2 ^=г?(г1 —1)2/2 ^3,°=г3(6^2 — 15зд+10) ^3,1=——4) ^1,°=-—4) 0,2 3 2 1,1 3 2,° 3 2 =z()z:L/2

ВП2 : ^1 = (0,1/2), ^=(1/2, 0), ^=(1/2,1/2)

^1 = 162221 г2/$11 Qз=—'&z°z{z1/f°l

Таблица 3

ИПЭ р=1, Ктах=4, к' =1, р' = 7

ВП1 : & = (0,1), 6=(1, 0),6 =(0,0)

°,° ¥>1 =—) °,° 4т> ! \ °,° = —24Р°(2°)

°,1 V 1 = 24(22-1)^1(^2) ^2,1=2422Р1(21) 0,1 =2422 Р1 (2°)

1,° Ч>1 = 2124 Р1 (22) ^2,° = ^4(^1 — 1)^1(21) 1,0 =2421Р1 (2°)

= — 24(22 — 1)2(422 —5)/2 ^2,2 = —(1/2)2422(421—5) 0,2 = — 2422(420 — 5)/2

1,1 = —21224(22 —1)(422 —5) ^1,1=—24 22(21 —1)(421—5) 1,1 =—242122(420—5)

2,° = —22224(422 —5)/2 ^2,° = —24 (21 —1)2(421—5)/2 2,° = — 242?(42о — 5)/2

0,3 = 24(22 —1)3/6 °,3 4 3 =2]_ 23 /6 °,3 =2423/6

1,2 = 2124(22 —1)2/2 ^1,2 = 2421(2l — 1)/2 1,2 =242122/2

2,1 ^ = 2224(22 —1)/2 ^2,1=2422(21 —1)2/2 2,1 = 242:122/2

3,° V 1 =2324/6 <^° = 24(21 —1)3/6 3,° =2423/6

ВП2 : = (0,1/2), ^=(1/2, 0), Ц3=(\/2,1/2)

$1 = 821 22 22 (1222 — 721 — 82122 — 822 +822 ) //и

Я2 = — 822 2222 (82!+82122 — 821+722 — 1222 )//22

Зз = 422 ^(^Ц — 1722+5 —1721+3221 г^^Ц)/^

ВП3 : С4=(1/4,1/2), Сб=(1/2,1/4), Сб=(1/4,1/4)

^4 = 102422 22 22 (422 — 1) ^5 = 10242222 22(421 —1) 36= 1024222222(42° —1)

р°(2 ^) = (202°—7022+842^ — 35), Р^) = (1022—242^+15)

4

ЗЩОШНЩОШНЩ01Ю

I ^ ^

2

соо

СОО ООО

I I

2

2

2

^ ^ ^

ООО ООО ООО ^ ^ ^

^ ^ ^ ^

ООО ООО ООО ООО ООО ^

^ ^ ^ ^ ^ ООО

^ « « « ^

ю ю со со со оо

оюою^юою^юс^ю

о

£

оо

X й 2

е ю (

О —1 <мсч

со со ОО

2) ^

^ 7 1 7

ГО (

К 2 00^ О0( 2

= ( 1 1 (

н

2

оооооо^ооооо^оос^оо

^ Со

I

2

I I

К! ООс^ 1

2 ООс^ 2 2

( ( ( (

ООс^ ООс^ ООс^ ООс^

1 1

+

+

Я

со^

со^

2

со^

05 05

со^

1

2 1

( ,—1

2

00^

О0(

ся ООс^ ООО

2

О0(

оо о 2

со

4+

оос^ ООО

I ^ -2

СЧ 4—( ^ ООС^ ООс^

с^ ^ ^

Ч—( ^ сч^ч О0(М ^ ^

^ ООО ООО ^

2

со

ООО

оо

2

со

++

о о

2

1 1 1

1 1

( \

2 2 2 т-Н

( ( 2 ( 2

ООс^ ООс^ ООс^ ООс^ О

О0( О0с(

00^ 00^ 00^ 00^

00^ 00^ 1 '

2 2 ^

со со оо оо

о 1= о 2

П

а. а. а. а. а. В

со С0

ся

ся

1 1 1

1 1 1

( ( )(

1 1 1

1 1 1

2 2 2

( ( )(

1 1 1

1 о 1 о 1

^

2 2 2

( ( (

^

^о с^ о о

^

оТ оТ с9

Ю

ю ю 2

( ( (

ся ^

с? с?

ся

ся о

1

о о 1 ся

)( )( )(

1

ся 1 ся

^

2 2 2

)( )( )(

1

о о 1

2 2 2

( ( (

СЯ 1—1 СЯ^н

N0 СЯ о о

^

оТ оТ

ю

Ю Ю 2

( ( (

=1 со

С? с? с?

ИПЭ р=1,

x=5, к'=2, р'=9

Таблица 5

ВП1 : & = (0,1), &=(1, 0), £з=(0, 0)

=z|To (Z2)

<^p=z0ro (¿1)

—0T0~

<^o' = 0,1 <^o' = 1'0 <^з' = 0'2 <^з' = 1'1 <^з' = 2'0 <^з' = 0'з <^з' = 1'2 <^з' = 2'1 <^з' =

^з'0=-^о^з(5^о-6)/6 0'4

<^з' = 1'з <^з' = 2'2 <^з' =

з'1 <^з' =

4'0 <^4' =

0'1 = -Z2°(Z2-1)T1(Z2)

1'0 = -Z1Z°T1(z2)

0'2

^0'2 = ^2(^2-1)2 t2(z2)/2 ^1'1=Z1 ^2(Z2-1)T2(Z2) ^1'0=z12zlT2(z2)/2 ^1'з=-г2(^2-1)з(5^2-6)/6 ^1'2=-Z1z1(z2-1)2(5z2-6)/2 ^2'1 = ^2 (^2 1)(5^2 6)/2 ^з'0=-^1з^2(5^2-6)/6 ^1'4=^2(^2-1)4/24

^'2=44(.2-1)2/4 ' 1 ^2 (^2 1)/6 ^1'0=z4z2/24

^2

0'1 = -^2^2T1(z1)

1 '0=-z2(z1-1)T1(z1)

^2

^2 ' 2=z2Z2T2(ZI )/2

' 1=Z2Z2(Z1-1)T2(Z1)

^2'0 = Z2(Z1-1)2T2(Z1 )/2

^0 ' з = -(5^1 -6)/6 1

2 ^2

1 ' 2=-z2z2(zi-1)(5zi —6)/2

2' 1=—z2z2(z1 — 1)2(5Z1—6)/2

^2

з'0=—Z2(Z1 —1)з(5^1—6)/6

^2'4=z24/24 = z0zQ(Z1-1)/6

=z2z22(z1 —1)2/4

= Z2Z2(Z1 — 1)з/6

0=^2(^1 —1)4/24

1 з

2 2

з 1

^2 '

: —Z2Z1T1(Z0) :Z2Z1Z2T2(Z0)

:4z2T2(z0)/2

:-z2zf(5z0-6)/6 :-z2z1z2(5z0-6)/2

:-z2Z2Z2(5^-6)/2

:44/24

:z02z4/24

ВП2 : ^ = (0,1/2), ^ = (1/2, 0), rço=(1/2,1/2), ^=(0,1/3), ^=(0, 2/3), ^6=(1/3, 0), ^r=(2/3, 0), ^8=(1/3, 2/3), rçg=(2/3,1/3)

^1=256^1 ^2 ((3^1 ^2 5^2 ^2+^2)^*11 4^1 (^2 2-0 )/12) 1

д4=(729/16)444(3зд-1)//?1

до = (729/16)444(3^-1)//121

^2 = 256^Q ((3^1 ^2 5 22 )/*22+4^2 (24 ^0 )^21 )/^22

дб = (729/16)2з2з22(320-1)//222 з з 2

Qo=

Q8 =

Qr=(729/16)z0z0z1 (3^1 2

128Z0 zq ((7z2 - 2z0 - Z1Z2 ) /01+2Z0 (Z1 - Z2 ) /02 ) //01

2„3.,3^ iw*2 Qg = (729/64)zgz0z1(3z1

: (729/64) z2Z3Z3(3Z2-1)//2

-1)//222

-1)//021

Q10 = 19683z3z3z1

ВП3 : Сю=(1/3,1/3)

T0(zj ) = (70zj-315zo+540zj-420zj+126)

ïi(zj ) = (35zO-120z2+140z2-56),Ï2(zi ) = (15z2-35Zj +21)

Замечание. При к' = 1 на равномерных сетках можно использовать базис с непрерывной первой производной, состоящий из приведённых ИПЭ (z') и Qs(z') при /ох = Уц = /22 = 1. При этом производные таких полиномов по нормали к границе в общем случае не удовлетворяют условиям (18).

Литература

1. Symbolic-Numerical Solution of Boundary-Value Problems with Self-Adjoint Second-Order Differential Equation using the Finite Element Method with Interpolation Hermite Polynomials / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky et al. // Lecture Notes in Computer Science. — 2014. — Vol. 8660. — Pp. 138-154.

2. Gusev A. A, Hai L. L., Chuluunbaatar O, Vinitsky S. I. Program KANTBP 4M for Solving Boundary-Value Problems for Systems of Ordinary Differential Equations of the Second Order. — http://wwwinfo.jinr.ru/programs/jinrlib/.

3. Сьярле Ф. Метод конечных элементов для эллиптических задач. — Москва: Мир, 1980.

4. Ramdas Ram-Mohan L. Finite Element and Boundary Element Applications in Quantum Mechanics. — New York: Oxford Univ. Press, 2002.

5. Algorithms for Solving the Parametric Self-Adjoint 2D Elliptic Boundary-Value Problem Using High-Accuracy Finite Element Method / A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, A. Gozdz // Вестник Российского университета дружбы народов. Серия: Математика, Информатика, Физика. — 2017. — Т. 25, № 1. — С. 36-55.

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

7. Шайдуров В. В. Многосеточные методы конечных элементов. — Москва: Наука, 1989.

8. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. — Москва: Стройиздат, 1982.

9. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — Москва: Наука, 1987.

10. Becker E. B., Carey G. F., Tinsley Oden J. Finite Elements. An Introduction. Vol. I. — New Jersey: Englewood Cliffs, Prentice Hall, 1981.

11. Стренг Г., Фикс Г. Теория метода конечных элементов. — Москва: Мир, 1977.

UDC 517.958, 530.145.6, 519.632.4 DOI: 10.22363/2312-9735-2017-25-3-217-233

High-Accuracy Finite Element Method for Solving Boundary-Value Problems for Elliptic Partial Differential

Equations

A. A. Gusev

Joint Institute for Nuclear Research 6 Joliot-Curie St., Dubna, Moscow region, 141980, Russian Federation

A new computational scheme of the finite element method of a high order of accuracy for solving boundary value problems for an elliptic partial differential equation that preserves the continuity of the derivatives of the approximate solution in a bounded domain of a multidimensional Euclidean space is proposed. A piecewise continuous basis of the finite element method is generated using interpolation Hermite polynomials of several variables and ensures the continuity of not only the approximate solution but also its derivatives up to a given order on the boundaries of finite elements, depending on the smoothness of the variable coefficients of the equation and the boundary of the domain. The efficiency and accuracy order of the computational scheme, algorithm and program are demonstrated by the example of an exactly solvable boundary-value problem for a triangular membrane depending on the number of finite elements of the partition of the domain and the dimension of the eigenvector of the algebraic problem. It was shown that, in order to achieve a given accuracy of the approximate solution, for schemes of the finite element method with Hermite interpolation polynomials the dimension of the eigenvector is approximately two times smaller than for schemes with Lagrange interpolation polynomials that preserve on the boundaries of finite elements only the continuity of the approximate solution. The high-accuracy computational scheme of the finite element method is oriented to calculations of the spectral and optical characteristics of quantum-mechanical systems.

Key words and phrases: elliptic partial differential equations, boundary-value problem, finite element method, interpolation polynomials

References

1. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, et al., Symbolic-Numerical Solution of Boundary-Value Problems with Self-Adjoint Second-Order Differential Equation using the Finite Element Method with Interpolation Hermite Polynomials, Lecture Notes in Computer Science 8660 (2014) 138-154.

2. A. A. Gusev, L. L. Hai, O. Chuluunbaatar, S. I. Vinitsky, Program KANTBP 4M for Solving Boundary-Value Problems for Systems of Ordinary Differential Equations of the Second Order.

URL http://wwwinfo.jinr.ru/programs/jinrlib/

3. P. Ciarlet, The Finite Element Method for Elliptic Problems, North-holland Publ. Comp, Amsterdam, 1978.

4. L. Ramdas Ram-Mohan, Finite Element and Boundary Element Applications in Quantum Mechanics, Oxford Univ. Press, New York, 2002.

5. A. A. Gusev, O. Chuluunbaatar, S. I. Vinitsky, V. L. Derbov, A. Gozdz, Algorithms for Solving the Parametric Self-Adjoint 2D Elliptic Boundary-Value Problem Using High-Accuracy Finite Element Method, RUDN Journal of Mathematics, Information Sciences and Physics 25 (1) (2017) 36-55.

6. O. A. Ladyzhenskaya, The Boundary Value Problems of Mathematical Physics, Vol. 49, Springer, Berlin, 1985.

7. V. V. Shaidurov, Multigrid Methods for Finite Elements, Springer, Berlin, 1995.

8. K. J. Bathe, Finite Element Procedures in Engineering Analysis, Englewood Cliffs, Prentice Hall, New York, 1982.

9. N. S. Bakhvalov, N. P. Zhidkov, G. M. Kobelkov, Numerical Methods, Nauka, Moscow, 1987, in Russian.

10. E. B. Becker, G. F. Carey, J. Tinsley Oden, Finite elements. An introduction. Vol. I, Englewood Cliffs, Prentice Hall, New Jersey, 1981.

11. G. Strang, G. J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, New York, 1973.

© Гусев А. А., 2017