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

Abstract


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.


1. Введение В работах [1, 2] разработаны символьно-численные алгоритмы и программы решения краевых задач для системы обыкновенных дифференциальных уравнений второго порядка методом конечных элементов (МКЭ) высокого порядка точности с применением интерполяционных полиномов Эрмита (ИПЭ). Такая реализация МКЭ в отличие от традиционной с использованием интерполяционных полиномов Лагранжа (ИПЛ) обеспечивает непрерывность производных до заданного порядка приближенного решения не только на конечных элементах сетки, но и на границах конечных элементов [3], т.е. обеспечивает законы сохранения тока носителей заряда в квантово-размерных полупроводниковых системах или тока вероятности в квантовомеханических задачах рассеяния [4]. Вычислительная схема МКЭ решения с высокой точностью краевых задач на собственные значения для эллиптического уравнения в частных производных в ограниченной области двумерного пространства с применением ИПЛ дана в работе [5]. Обобщение вычислительной схемы МКЭ для решения эллиптических краевых задач с применением ИПЭ составляет цель данной работы. В настоящей работе представлена вычислительная схема МКЭ решения с высокой точностью эллиптических краевых задач в ограниченной области многомерного евклидова пространства, заданного в виде многогранника. Кусочно-непрерывный базис МКЭ генерируется с помощью ИПЭ нескольких переменных и обеспечивает непрерывность не только приближенного решения, но и его производных до заданного порядка в зависимости от гладкости переменных коэффициентов уравнения и границы области. Эффективность вычислительной схемы, алгоритма и программы демонстрируется на примере эталонной точно-решаемой краевой задачи на собственные значения для треугольной мембраны. Структура работы следующая. В разделе 2 дана постановка задачи. В разделе 3 дана формулировка вычислительной схемы МКЭ с применением ИПЭ нескольких переменных. В разделе 4 представлены результаты расчётов эталонной краевой задачи, демонстрирующие эффективность вычислительной схемы и порядок точности по числу конечных элементов разбиения области и числу кусочно-полиномиальных функций при различных ИПЭ двух переменных. В заключительном разделе даны перспективы развития и применения предложенных вычислительных схем. 2. Постановка задачи Рассмотрим самосопряжённую краевую задачу для эллиптического дифференциального уравнения второго порядка: Для коэффициентов главной части (1) выполняется условие равномерной эллиптичности в ограниченной области =(1, . . . , ) Ω евклидова пространства , т.е. существуют константы > 0, > 0 такие, что Также предполагается, что 0() > 0, () = () и () - функции, непрерывные вместе со своими обобщёнными производными до заданного порядка в области, ∈ Ω¯ = Ω ∪ ∂Ω с кусочно-полиномиальной границей = ∂Ω, обеспечивают существование нетривиальных решений, подчинённых граничным условиям первого (I), второго (II) или третьего (III) рода [6, 7]: где () = (()) - известная функция на границе области () , ∂Φ()/∂ - производная по направлению конормали, ˆ - внешняя нормаль к границе области ∑︀ = ∂Ω, ˆ - единичный вектор вектора = ˆ, (ˆ, ˆ) - скалярное произведение в . Для задачи дискретного спектра собственные функции Φ() из пространства Соболева 1(Ω), Φ() ∈ 1(Ω), соответствующие собственным значениям энергии : 1 � 2 � . . . � � . . . удовлетворяют условиям нормировки и ортогональности Решение МКЭ краевых задач (1)-(3) сводится к нахождению стационарных точек вариационного функционала [7, 8] где Π(Φ, , ) - симметричный квадратичный функционал вида (5) Поверхностный интеграл в (4) вычисляется с учётом (2). 3. Вычислительная схема МКЭ В МКЭ область Ω = Ω() = ∆ , заданная в виде многогранника, покрывается конечными элементами, в данном случае симплексами ∆ с + 1 вершинами ˆ = (ˆ1, ˆ2, . . . , ˆ) при = 0, . . . , . Каждое ребро симплекса ∆ разбиваем на равных частей и проводим семейства параллельных гиперплоскостей (, ), нумеруя каждую целым числом от = 0, . . . , , начиная от соответствующей грани, например, как показано при = 2 на рис. 1 (см. также [9, С. 220]). Уравнение гиперплоскости (, ) : (; ) - / = 0, где (; ) - линейная функция от . Узловые точки пересечения гиперплоскостей нумеруем наборами целых чисел [0, . . . , ], ): 0, 0 + . . . + = , где , = 0, 1, . . . , - номера гиперплоскостей, параллельных грани симплекса, не содержащей -вершину ˆ = (ˆ1, . . . ˆ). Координаты = (1, . . . , ) узловой точки ∈ ∆ вычисляются по формуле (1, . . . , ) = (ˆ01, . . . , ˆ0)0/ + (ˆ11, . . . , ˆ1)1/ + . . . + (ˆ1, . . . , ˆ)/ (6) через координаты вершин ˆ = (ˆ1, . . . , ˆ). Тогда ИПЛ () равные единице в точке с координатами = (1, . . . , ), характеризуемой числами [0, 1, . . . , ], и нулю в остальных точках ′ , т.е. (′ ) = ′ имеют вид Рис. 1. (a) Нумерация узлов , = 1, . . . , ( + 1)( + 2)/2, наборами чисел [0, 1, 2] стандартного треугольного элемента Δ для схемы с ИПЛ пятого порядка ′ = = 5 при = 2; линии (пять пересекающихся прямых) - нули ИПЛ 14(′) из (13), равного единице в точке, пронумерованной тройкой чисел [0, 1, 2] = [2, 2, 1]; (b) изолинии ИПЛ 14(′) Отметим, что построение ИПЭ () с фиксированными значениями функций {(′ )} и производных {∂∙()|= ′ } в узлах ′ уже при = 2 приводит к громоздким выражениям, затратным для использования в МКЭ на неравномерных сетках. Экономичная реализация, принятая в МКЭ, следующая: § расчёты проводятся в локальных координатах ′, в которых координаты вершин симплекса ∆ следующие: ˆ′ = (ˆ′ 1, . . . , ˆ′ ), ˆ′ = , § ИПЭ в исходных координатах ищется в виде линейных комбинаций полиномов в локальных координатах ′, при этом переход в исходные координаты выполняется только на этапе численного решения конкретной краевой задачи (1)-(4), § вычисление интегралов МКЭ выполняется в локальных координатах. Построим ИПЭ на произвольном -мерном симплексе ∆ с + 1 вершинами ˆ = (ˆ1, ˆ2, . . . , ˆ), = 0, . . . , . Для этого введём локальную систему координат ′ = (1′ , 2′ , . . . , ′ ) , в которой координаты вершин симплекса следующие: ˆ′ = (ˆ′ = , = 1, . . . , ). Связь между координатами даётся формулой: , (8) =1 или в матричном виде (9) Обратное преобразование имеет вид . (10) Отсюда следует связь между операторами дифференцирования Формула (11) применяется для вычисления ИПЭ (′)={ˇ(′), (′)} из (21), удовлетворяющих условиям (14), (18), (19) из следующего раздела, с фиксированными производными до заданного порядка в узлах ′ . При этом производные по нормали к границе элемента в исходной системе координат в общем случае не являются таковыми в локальных координатах ′. При построении ИПЭ в локальных координатах ′ требуется пересчёт фиксированных производных в узлах ′ элемента ∆ в узлы ′ ′ элемента ∆, используя матрицы ˆ-1, заданные громоздкими выражениями. Поэтому требуемый пересчёт производных выполняется, используя соотношения (8)-(11), для каждого конечного элемента ∆ на этапе формирования базиса ИПЭ {¯′ (′)} на конечном элементе ∆, который выполняется численно по аналитическим формулам, представленным в следующем разделе. Интегралы, входящие в вариационный функционал (4) на области Ω() = , выражаем через интегралы, вычисленные на элементе ∆ , пересчитанные в локальные координаты ′ на элементе ∆, (12) где = det ˆ>0 определитель матрицы 1. Интерполяционные полиномы Лагранжа В локальных координатах ИПЛ (′), равные единице в узловой точке ′ , характеризуемой числами [0, 1, . . . , ], и нулю в остальных узловых точках ′ ′ , т.е., (′ ′ ) = ′ определяются по формуле (7) при (0; ′)=1 1′ . . . ′ , (; ′)=′, . (13) Приравнивание числителей в (13) нулю даёт семейства уравнений прямых, направленных «горизонтально», «вертикально» и «наклонно» в локальной системе координат элемента ∆, которая связана аффинным преобразованием (8) с семействами «наклонно» расположенных прямых элемента ∆ . На рис. 1 дан пример, иллюстрирующий построение ИПЛ при = 2, , ′ = 1, . . . , ( + 1)( + 2)/2, = 5 на элементе ∆ в виде прямоугольного треугольника с вершинами ˆ0′ = (ˆ0′ 1, ˆ0′ 2) = (0, 0), ˆ1′ = (ˆ1′ 1, ˆ1′ 2) = (1, 0), ˆ2′ = (ˆ2′ 1, ˆ2′ 2) = (0, 1). Кусочно-полиномиальные функции ¯(), формирующие конечноэлементный базис {¯()} , которые строятся путём сшивки ИПЛ () из (7), получаемых из 1. с помощью преобразования (10), на конечных элементах ∆ : () = {(), ∈ ∆ ; 0, ̸∈ ∆ } , являются непрерывными, но их производные терпят разрывы на границах элементов ∆ . 2. Интерполяционные полиномы Эрмита Построим ИПЭ порядка ′, сшивкой которых можно получить кусочнополиномиальные функции (25), сохраняющие непрерывность производных до заданного порядка ′. 1. Вспомогательные полиномы (ВП1) Для построения ИПЭ в локальных координатах ′ введём набор вспомогательных полиномов (ВП1) Здесь в узловых точках ′ , определённых согласно (6), в отличие от ИПЛ заданы значения не только функций, но и их производных до порядка max 1. ВП1 даются выражениями где коэффициенты 1... ,1... вычисляются из рекуррентных соотношений, полученных в результате подстановки (15) в условия (14), При > 1 и max > 1, число max′ ИПЭ степени ′ и кратности узлов max меньше числа полиномов 1′ , формирующих базис в пространстве полиномов степени ′ (например, ИПЛ из (13)), т.е. полиномы, удовлетворяющие (14), определены неоднозначно. 3.2.2. Вспомогательные полиномы (ВП2 и ВП3) Для однозначного определения полиномиального базиса введём = 1′ max′ вспомогательных полиномов () двух типов: ВП2 и ВП3, линейно-независимых от ВП1 из (15) и удовлетворяющих условиям в узловых точках ′ ′ ВП1: , (18) Отметим, что для обеспечения непрерывности производных часть полиномов, называемых ВП2, должны удовлетворять условию где ′ ′ = (′ ′1, . . . , ′ ′) - выбранные точки, лежащие на гранях всевозможных размерностей (от 1 до 1) -мерного симплекса ∆ и не совпадающие с узловыми точками ИПЭ ′ , в которых выполняется (14), ∂/∂() - производная по направлению вектора , направленного по нормали к соответствующей -й грани -мерного симплекса ∆ в точке ′ в исходной системе координат, который пересчитывается в точку ′ ′ грани симплекса ∆ в локальной системе координат, используя соотношения (8)-(11), например, при = 2 см. (23). Вычисляя число 1() независимых параметров, требуемых для обеспечения непрерывности производных до порядка , определяем его максимальное значение ′, которое можно получить для схем с заданными и max, и, соответственно, дополнительных условий (18). Остаётся 2= 1(′) независимых параметров и, соответственно, добавляется 2 дополнительных условий, необходимых для однозначного определения полиномов, называемых ВП3, (19) где ′′ = (′′1, . . . , ′′) ∆ выбранные точки, принадлежащие симплексу без границы, но не совпадающие с узловыми точками ВП1 ′ . Вспомогательные полиномы ВП2 даются выражением (20) где = 1, если точка , в которой заданы дополнительные условия (18), лежит на соответствующей грани симплекса ∆, т.е. (, ) = 0, и = ′, если (, ) = 0. Вспомогательные полиномы ВП3 даются выражением (20) при = ′. Коэффициенты 1,..., ; определяются из системы линейных уравнений, полученных в результате подстановки (20) в условия (17)-(19). В результате получаем требуемый набор базисных ИПЭ , (21) составленных из полиномов (′) типа ВП2 и ВП3 и полиномов ˇ(′) типа ВП1: 3.2.3. ИПЭ при = 2 При = 2 степень ′ полинома по тангенциальной переменной на границе треугольника ∆ совпадает со степенью полинома двух переменных, и для его однозначного определения требуется ′ + 1 параметр. Производная порядка ′ по нормальной переменной к границе будет полином степени ′ ′, и для её однозначного определения потребуется ′ ′ +1 параметров. Однако, она определена только ′ ′(+1) параметрами: смешанными производными фиксированного порядка ′ по нормальной переменной и порядка от 0 до max ′ 1 по тангенциальной переменной, т.е. Таким образом, с учётом того, что у треугольника три стороны, для однозначного определения производных до порядка ′ на границе необходимо 1(′) = 3 + . . . + 3′ = 3′(′ + 1)/2 параметров и, соответственно, дополнительных условий (18). Например, для = 1 и max = 4 имеется = 6 дополнительных условий для определения ВП2 и ВП3. Степень ′ = 7 полинома по тангенциальной переменной на границе треугольника совпадает со степенью полинома двух переменных, и для его однозначного определения требуется ′ + 1 = 8 параметров. Производная первого порядка ′ = 1 по нормальной переменной к границе будет полином степени ′ ′ = 6, и для её однозначного определения потребуется ′ ′ +1 = 7 параметров. Однако она определена только ′ ′( + 1) = 6 параметрами: смешанными производными ∂/∂, ∂2/(∂∂), ∂3/(∂∂2), заданными в двух вершинах. Недостающий параметр можно определить, задавая производную по направлению, непараллельном границе треугольника в одной из точек на его стороне (например, посередине). Таким образом, для = 1 и max = 4 можно построить ИПЭ с фиксированными значениями первой производной на границе треугольника, при этом остаётся 6 3 = 3 свободных параметров. Производная второго порядка ′ = 2 по нормальной переменной к границе - полином степени ′ - ′ = 5, и для её однозначного определения потребуется ′ ′ + 1 = 6 параметров. Однако она определена только ′ ′( + 1) = 4 параметрами: смешанными производными ∂2/∂2 и ∂3/(∂2∂), заданными в двух вершинах треугольника. Таким образом, для однозначного определения второй производной потребуется 6 параметров. Это означает, что данным алгоритмом для = 1 и max = 4 нельзя построить схему МКЭ с непрерывной второй производной. В этом случае следует использовать схему с max > 4, например, обозначенную [152] из табл. 1 и рис. 2. Тогда оставшиеся три свободных параметра используются для построения ВП3. Заметим, что возможно построение схем, обеспечивающих непрерывность вторых производных на некоторых границах конечных элементов. Этот случай в данной работе не рассматривается. При = 2 производные ∂/∂ по направлению , перпендикулярном к соответствующей грани = 0, 1, 2 в исходной системе координат, даются через частные производные ∂/∂′ , = 1, 2 в локальной системе координат ∆, используя (8)-(11), выражениями (23) где = (ˆ0, ˆ1, ˆ2) - функции от координат вершин ˆ0, ˆ1, ˆ2 треугольника ∆ в исходной системе координат Схематическое изображение реализации условий (14), (17), (18) и (19), с помощью которых построены базисные ИПЭ, показано при = 2 на рис. 2. Характеристики полиномиального базиса из ИПЭ на элементе ∆ при = 2 приведены в табл. 1. Выражения для ИПЭ (21) до ′ = 9 порядка при = 2 представлены в Приложении. 3.2.4. Кусочно-полиномиальные функции Кусочно-полиномиальные функции () с непрерывными производными до порядка ′ строятся путём сшивки полиномов () = {ˇ(), ()} из (21), получаемых алгоритмом из Приложения, на конечных элементах ∆ ∈ Ω() = ⋃︀ ∆ : , (25) где знак «-» может появиться только для ВП2, когда нужно сшить нормальные производные нечётного порядка. Рис. 2. Схематическое изображение условий на элементе Δ для построения базиса ИПЭ [max′]: [131], [141], [231], [152]. Квадраты: точки ′ , в которых фиксируются значения функций и их производных согласно условиям (14), (17); сплошные (прерывистые) стрелки: с началом в точках ′ , в которых фиксируются значения первой (второй) производной по направлению нормали в исходных координатах, согласно условию (18), соответственно; кружки: точки ′ , в которых фиксируются значения функций согласно условию (19) Характеристики базиса ИПЭ (21) при = 2 Таблица 1 [max′] [131] [141] [231] [152] ′ max( + 1) - 1 5 7 8 9 max′ ( + 1)( + 2)max(max + 1)/4 18 30 36 45 1′ (′ + 1)(′ + 2)/2 21 36 45 55 ( + 1)max(max - 1)/4 3 6 9 10 1(1) 3 3 3 6 3 1(2) 9 9 9 18 9 (ВП1) max′ 18 30 36 45 (ВП2) 1(′) 3 3 6 9 (ВП3) - 1(′) 0 3 3 1 Ограничение на порядок производных ′ : 3′(′ + 1)/2 � . Разложение искомого решения Φ() по базису кусочно-полиномиальных функций ′ (), Φ () = ∑︀ ′ ()Φ и подстановка его в вариационный функционал (4) приводит к обобщённой алгебраической задаче на собственные значения, ( - )Φ = 0 , которая решается стандартным методом (см., например, [8]). Элементы симметричных матриц жёсткости и масс содержат интегралы типа (4), которые вычисляются на элементах в области ∆ ∈ Ω() = ⋃︀ ∆ , пересчитанные в локальные координаты на элементе ∆. Теоретические оценки погрешностей между точным , Φ() ∈ 2(Ω) и приближенным , Φ () ∈ ′+1 1(Ω) решениями следующие [10, 11]: , (26) где - максимальный размер конечного элемента ∆ , ′ - порядок схемы МКЭ, - номер собственного значения, 1 и 2 - коэффициенты, не зависящие от . 4. Результаты и обсуждение В качестве примера рассмотрим решение задачи дискретного спектра (1)-(3) при = 2, 0() = () = 1, и () = 0 в области Ω() = ∪=1∆ в виде равностороннего треугольника со стороной 4/3 с граничными условиями второго типа (II), который разбит на = 2 равносторонних треугольников ∆ со стороной = 4/3. Собственные значения этой задачи с вырожденным спектром - целые числа = 0, 1, 1, 3, 4, 4, 7, 7, . . .. 4 На рис. 3 представлены конечноэлементная сетка с ИПЛ восьмого порядка и профиль четвёртой собственной функции Φ(). Рис. 3. (a) Сетка на области Ω() = ⋃︀ Δ треугольной мембраны, составленная из треугольных элементов Δ ; 0. профили четвёртой собственной функции Φ() = 3, полученной на ИПЛ порядка ′ = = 8 На рис. 4 показаны погрешности ∆ = - собственного значения () в зависимости от числа (число элементов 2) и от длины вектора Φ алгебраической задачи на собственные значения для схем МКЭ от пятого до девятого порядка точности, использующие ИПЛ, отмеченные метками [max′] = [510], . . . , [910], и использующие ИПЭ, отмеченные метками [131], [141], [231] и [152] из табл. 1, сохраняющие непрерывность первой и второй производных приближённого решения соответственно. Как видно из рис. 4, погрешности собственного значения ∆() схем МКЭ одного порядка ′ = max( + 1) 1 примерно одинаковые и соответствуют теоретическим оценкам (26), но в схемах МКЭ с ИПЭ, сохраняющих непрерывность первой и второй производных приближённого решения, используются матрицы меньшей размерности, соответствующие длине вектора , в 1,5-2 раза меньшей, чем для схем с ИПЛ, сохраняющих на границах конечных элементов только непрерывность приближённого решения. Вычисления проводились на компьютере 2 x Xeon 3.2 GHz, 4 GB RAM, используя Intel Fortran 77 с четверной точностью real*16, с 32 значащими цифрами. Время счёта рассмотренных примеров не превышало 3 минут. Рис. 4. Зависимости погрешности Δ собственного значения от числа элементов и от длины вектора 5. Заключение Предложена вычислительная схема МКЭ высокого порядка точности решения задачи на собственные значения для эллиптического уравнения в частных производных в ограниченной области многомерного евклидова пространства, обеспечивающая непрерывность не только приближенного решения, но и его производных до заданного порядка. На примере точно-решаемой краевой задач для треугольной мембраны показано, что для достижения заданной точности приближённого решения, для схем с МКЭ с ИПЭ, обеспечивающих непрерывность первой и второй производных приближённого решения, используются матрицы меньшей размерности, соответствующие длине вектора , в 1.5-2 раза меньшей, чем для схем с ИПЛ, сохраняющих на границах конечных элементов только непрерывность приближённого решения. Отметим, что явное выражение для ИПЭ нескольких переменных, в отличие от ИПЛ, достаточно сложное, в локальных координатах ′ требуется пересчёт фиксированных производных в узлах ′ элемента ∆ в узлы ′ ′ элемента ∆, используя матрицы ˆ-1. Поэтому требуемый пересчёт производных выполняется для каждого конечного элемента ∆ на этапе формирования базиса ИПЭ {¯′ (′)} на конечном элементе ∆, который выполняется численно по аналитическим формулам, аналогичным формулам при = 2 из раздела 3.4 и Приложения. Вычислительные схемы МКЭ ориентированы на расчёты спектральных и оптических характеристик квантовых точек и других квантовомеханических систем. Реализация МКЭ c ИПЭ в пространстве ): 2 и других областях, отличных от многогранника, будет дана в последующих работах. 6. Приложение Выражения для ИПЭ в исходных координатах, удовлетворяющие условиям (14), (18), (19), так, что для каждого полинома в правых частях этих условий единица появляется только один раз, громоздкие. Их вычисление удобно проводить по алгоритму: 1) по формулам (15) и (20) вычисляются ВП1 (′), ВП2 и ВП3 (′); 1. используя формулу (22), вычисляем базис ИПЭ, состоящий из полиномов ˇ(′) и (′), в локальной системе координат; 2. ВП1 () пересчитываются по формулам (11); 3. выполняется преобразование (10). В табл. 2-5 приведены ВП1 (′), ВП2 и ВП3 (′), а соответствующие коэффициенты ;; вычисляются по формулам (22). Обозначения: , , - координаты узлов в которых правая часть (14), (18) или (19) равна единице, 0 = 1 1 2, определяются из формул (24), при этом аргументы у функций и штрихи у независимых переменных опущены. ИПЭ = 1, max = 3, ′ = 1, ′ = 5 (элемент Аргириса [3]) Таблица 2 ВП1 : 1=(0, 1), 2=(1, 0), 3=(0, 0) 0,0=3(62-152+10) 1 2 2 0,1=-3(2-1)(32-4) 1 2 1,0=-13(32-4) 1 2 0,2=3(2-1)2/2 1 2 1,1=13(2-1) 1 2 2,0=23/2 1 1 2 0,0=3(62-151+10) 2 1 1 0,1=-32(31-4) 2 1 1,0=-3(1-1)(31-4) 2 1 0,2=32/2 2 1 2 1,1=(1-1)32 2 1 2,0=3(1-1)2/2 2 1 0,0=3(62-150+10) 3 0 0 0,1=-32(30-4) 3 0 1,0=-31(30-4) 3 0 0,2=32/2 3 0 2 1,1=312 3 0 2,0=32/2 3 0 1 ВП2 : 1=(0, 1/2), 2=(1/2, 0), 3=(1/2, 1/2) 1=16212/11 0 2 2=16222/22 0 1 2 2 3=-801 2 /01 ИПЭ =1, max=4, ′=1, ′=7 Таблица 3 ВП1 : 1=(0, 1), 2=(1, 0), 3=(0, 0) 0,0=-40(3) 1 2 0,1=4(2-1)1(2) 1 2 1,0=141(2) 1 2 0,2=-4(2-1)2(42-5)/2 1 2 1,1=-14(2-1)(42-5) 1 2 2,0=-24(42-5)/2 1 1 2 0,3=4(2-1)3/6 1 2 1,2=14(2-1)2/2 1 2 2,1=24(2-1)/2 1 1 2 3,0=34/6 1 1 2 0,0=-40(1) 2 1 0,1=421(1) 2 1 1,0=4(1-1)1(1) 2 1 0,2=-(1/2)42(41-5) 2 1 2 1,1=-42(1-1)(41-5) 2 1 2,0=-4(1-1)2(41-5)/2 2 1 0,3=43/6 2 1 2 1,2=42(1-1)/2 2 1 2 2,1=42(1-1)2/2 2 1 3,0=4(1-1)3/6 2 1 0,0=-40(0) 3 0 0,1=421(0) 3 0 1,0=411(0) 3 0 0,2=-42(40-5)/2 3 0 2 1,1=-412(40-5) 3 0 2,0=-42(40-5)/2 3 0 1 0,3=43/6 3 0 2 1,2=412/2 3 0 2 2,1=422/2 3 0 1 3,0=43/6 3 0 1 ВП2 : 1=(0, 1/2), 2=(1/2, 0), 3=(1/2, 1/2) 1=8122(122-71-812-82+82)/11 2 0 1 2 2 2 2 2 2=-81 20 (81 +812-81+72-122 )/22 3=4220(122-172+5-171+3212+122)/01 1 2 2 1 ВП3 : 4=(1/4, 1/2), 5=(1/2, 1/4), 6=(1/4, 1/4) 4=1024222(42-1) 0 1 2 5=1024222(41-1) 0 1 2 6=1024222(40-1) 0 1 2 0( ) = (203-702+84 -35), 1( ) = (102-24 +15) ИПЭ =2, max=3, ′=1, ′=8 Таблица 4 ВП1 : 1=(0, 1), 2=(1/2, 1/2), 3=(1, 0)4=(0, 1/2), 5=(1/2, 0), 6=(0, 0) 0,0=3(22-1)30(2) 1 2 0,0=3(21-1)30(1) 3 1 0,0=3(20-1)30(0) 6 0 0,1=-3(2-1)1(2) 1 2 0,1=-321(1) 3 1 0,1=-321(0) 6 0 1,0=-131(2) 1 2 1,0=-3(1-1)1(1) 3 1 1,0=-311(0) 6 0 0,2=3(2-1)2(22-1)3/2 1 2 0,2=3(21-1)32/2 3 1 2 0,2=32(20-1)3/2 6 0 2 1,1=3(22-1)31(2-1) 1 2 1,1=32(1-1)(21-1)3 3 1 1,1=312(20-1)3 6 0 2,0=3(22-1)32/2 1 2 1 2,0=3(1-1)2(21-1)3/2 3 1 2,0=32(20-1)3/2 6 0 1 0,0=64332(0) 2 1 2 0,1=3233(22-1)(60+1) 2 1 2 1,0=3233(21-1)(60+1) 2 1 2 0,2=833(22-1)2 2 1 2 1,1=1633(21-1)(22-1) 2 1 2 2,0=833(21-1)2 2 1 2 0,0=64332(1) 4 0 2 0,1=3233(22-1)(61+1) 4 0 2 1,0=64313(61+1) 4 0 2 0,2=833(22-1)2 4 0 2 1,1=32313(22-1) 4 0 2 2,0=32323 4 0 1 2 0,0=64332(2) 5 0 1 0,1=64332(62+1) 5 0 1 1,0=3233(21-1)(62+1) 5 0 1 0,2=32332 5 0 1 2 1,1=32332(21-1) 5 0 1 2,0=833(21-1)2 5 0 1 ВП2 : 1=(0, 1/4), 2=(0, 3/4), 3=(1/4, 0), 4=(3/4, 0), 5=(1/4, 3/4), 6=(3/4, 1/4) 1= (512/9)212(20-1)(22-1)(40-1)/11 0 2 2=-(512/9)212(20-1)(22-1)(42-1)/11 0 2 3=-(512/9)222(20-1)(21-1)(40-1)/22 0 1 4=-(512/9)222(20-1)(21-1)(41-1)/22 0 1 5= (256/9)022(21-1)(22-1)(42-1)/01 1 2 6= (256/9)022(21-1)(22-1)(41-1)/01 1 2 ВП3 : 7=(1/4, 1/2), 8=(1/2, 1/4), 9=(1/4, 1/4) 0 1 2 0 1 2 0 1 2 0( ) = (482-1052+58), 1( ) = (2 -1)3(9 -10), 2( ) = (242-12012/ +4) 2 ИПЭ =1, max=5, ′=2, ′=9 Таблица 5 ВП1 : 1=(0, 1), 2=(1, 0), 3=(0, 0) 0,0=50(2) 1 2 0,1=-5(2-1)1(2) 1 2 1,0=-151(2) 1 2 0,2=5(2-1)22(2)/2 1 2 1,1=15(2-1)2(2) 1 2 2,0=252(2)/2 1 1 2 0,3=-5(2-1)3(52-6)/6 1 2 1,2=-15(2-1)2(52-6)/2 1 2 2,1=-25(2-1)(52-6)/2 1 1 2 3,0=-35(52-6)/6 1 1 2 0,4=5(2-1)4/24 1 2 1,3=15(2-1)3/6 1 2 2,2=25(2-1)2/4 1 1 2 3,1=35(2-1)/6 1 1 2 4,0=45/24 1 1 2 0,0=50(1) 2 1 0,1=-521(1) 2 1 1,0=-5(1-1)1(1) 2 1 0,2=522(1)/2 2 1 2 1,1=52(1-1)2(1) 2 1 2,0=5(1-1)22(1)/2 2 1 0,3=-53(51-6)/6 2 1 2 1,2=-52(1-1)(51-6)/2 2 1 2 2,1=-52(1-1)2(51-6)/2 2 1 3,0=-5(1-1)3(51-6)/6 2 1 0,4=54/24 2 1 2 1,3=53(1-1)/6 2 1 2 2,2=52(1-1)2/4 2 1 2 3,1=52(1-1)3/6 2 1 4,0=5(1-1)4/24 2 1 0,0=50(0) 3 0 0,1=-521(0) 3 0 1,0=-511(0) 3 0 0,2=522(0)/2 3 0 2 1,1=5122(0) 3 0 2,0=522(0)/2 3 0 1 0,3=-53(50-6)/6 3 0 2 1,2=-512(50-6)/2 3 0 2 2,1=-522(50-6)/2 3 0 1 3,0=-53(50-6)/6 3 0 1 0,4=54/24 3 0 2 1,3=513/6 3 0 2 2,2=522/4 3 0 1 2 3,1=532/6 3 0 1 4,0=54/24 3 0 1 ВП2 : 1=(0, 1/2), 2=(1/2, 0), 3=(1/2, 1/2), 4=(0, 1/3), 5=(0, 2/3), 6=(1/3, 0), 7=(2/3, 0), 8=(1/3, 2/3), 9=(2/3, 1/3) 1=256313((312-52-2+2)11-41(2-0)12)/ 2 0 2 1 2 11 4=(729/16)323(30-1)/ 2 5=(729/16)323(32-1)/ 2 0 1 2 11 0 1 2 11 2=256332((312-52-2+1)22+42(1-0)21)/ 2 0 1 2 1 22 6=(729/16)332(30-1)/ 2 7=(729/16)332(31-1)/ 2 0 1 2 22 0 1 2 22 3=128033((72-20-12)01+20(1-2)02)/ 2 1 2 0 01 8=(729/64)233(32-1)/ 2 9=(729/64)233(31-1)/ 2 0 1 2 01 0 1 2 01 ВП3 : 10=(1/3, 1/3) 10=19683333 0 1 2 0( ) = (704-3153+5402-420 +126) 1( ) = (353-1202+1402-56), 2( ) = (152-35 +21) 2 2 Замечание. При ′ = 1 на равномерных сетках можно использовать базис с непрерывной первой производной, состоящий из приведённых ИПЭ ˇ(′) и (′) при 01 = 11 = 22 = 1. При этом производные таких полиномов по нормали к границе в общем случае не удовлетворяют условиям (18).

A A Gusev

Joint Institute for Nuclear Research

Author for correspondence.
Email: gooseff@jinr.ru
6 Joliot-Curie St., Dubna, Moscow region, 141980, Russian Federation

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

  • 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.
  • 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.
  • P. Ciarlet, The Finite Element Method for Elliptic Problems, North-holland Publ. Comp, Amsterdam, 1978.
  • L. Ramdas Ram-Mohan, Finite Element and Boundary Element Applications in Quantum Mechanics, Oxford Univ. Press, New York, 2002.
  • A.A. Gusev, O. Chuluunbaatar, S.I. Vinitsky, V.L. Derbov, A.G´o´zd´z, 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.
  • O.A. Ladyzhenskaya, The Boundary Value Problems of Mathematical Physics, Vol. 49, Springer, Berlin, 1985.
  • V.V. Shaidurov, Multigrid Methods for Finite Elements, Springer, Berlin, 1995.
  • K.J. Bathe, Finite Element Procedures in Engineering Analysis, Englewood Cliffs, Prentice Hall, New York, 1982.
  • N.S. Bakhvalov, N.P. Zhidkov, G.M. Kobelkov, Numerical Methods, Nauka, Moscow, 1987, in Russian.
  • E.B. Becker, G.F. Carey, J. Tinsley Oden, Finite elements. An introduction. Vol. I, Englewood Cliffs, Prentice Hall, New Jersey, 1981.
  • G. Strang, G.J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, New York, 1973.

Views

Abstract - 1446

PDF (Russian) - 100


Copyright (c) 2017 Gusev A.A.

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.