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

Полный текст

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).

Об авторах

Александр Александрович Гусев

Объединённый институт ядерных исследований

Автор, ответственный за переписку.
Email: gooseff@jinr.ru

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

ул. Жолио-Кюри, д. 6, г. Дубна, Московская область, Россия, 141980

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