CC BY
78
УДК 517.958:530.145.6
Вычислительные схемы для решения задачи Штурма—Лиувилля методом конечных элементов с интерполяционными полиномами Эрмита
A. A. Гусев*, Л. Л. Хай
* Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московской обл., Россия, 141980
^ Белгородский государственный национальный исследовательский университет ул. Победы, д. 85, г. Белгород, Россия, 308015
Построены вычислительные схемы решения задачи Штурма-Лиувилля с однородными краевыми условиями первого, второго и третьего рода методом конечных элементов, сохраняющие в приближённых решениях свойства непрерывности производных искомых решений. Выведены рекуррентные соотношения для вычисления в аналитическом виде интерполяционных полиномов Эрмита с узлами произвольной кратности. Из интерполяционных полиномов Эрмита сконструированы базисные кусочно—полиномиальные функции на конечноэлементной сетке с переменным шагом, аппроксимирующие решение исходной задачи. Исходная задача Штурма-Лиувилля в базисе кусочно полиномиальных функций редуцируется к обобщённой алгебраической задаче на собственные значения с ленточными матрицами жёсткости и масс. Построены матрицы жёсткости и масс в виде сумм интегралов, содержащих заданные коэффициентные и потенциальные функции исходного самосопряжённого дифференциального уравнения и вычисленные интерполяционные полиномы Эрмита и их производные. Интегрирование выполняется с помощью гауссовых квадратур, а в специальных случаях, включающих кусочно-полиномиальные коэффициентные и потенциальные функции, в аналитическом виде. Эффективность и скорость сходимости предложенных вычислительных схем и разработанных алгоритмов и программ в среде Maple-Fortran доказана численным анализом тестовых расчётов точно решаемых задач Штурма-Лиувилля с непрерывными и кусочно-непрерывными потенциальными функциями.
Ключевые слова: задача Штурма—Лиувилля, вычислительная схема, метод конечных элементов, интерполяционные полиномы Эрмита.
1. Введение
Изучение математических моделей, описывающих процессы распространения мод в плавно-нерегулярных волноводных системах, туннелирования и каналиро-вания составных квантовых систем через многомерные барьеры, фотоионизации и фотоабсорции в молекулярных, атомных и квантово-размерных полупроводниковых системах, требует разработки высокоточных и экономичных алгоритмов и программ решения краевых задач [1—11].
В этом направлении в рамках вариационно-проекционных формулировок краевых задач и метода конечных элементов (МКЭ), использующих интерполяционные лагранжевы элементы [12-14], были разработаны и применены для решения ряда квантовомеханических задач вычислительные схемы, символьно-численные алгоритмы и программы [7,15, 16]. Соответствующая алгоритмическая и программная реализация МКЭ, использующая лагранжевы элементы, сохраняла только свойство непрерывности искомого решения при его аппроксимации на конечноэлементной сетке. Однако для решения данного класса задач, особенно описывающих процессы в квантоворазмерных полупроводниковых системах и плавно-нерегулярных волноводах с кусочно-непрерывными коэффициентными
Статья поступила в редакцию 18 сентября 2014 г. Авторы благодарят С. И. Виницкого, В. П. Гердта, В. Л. Дербова, Л. А. Севастьянова и О. Чулуунбаатара за сотрудничество и поддержку. Работа поддержана грантами РФФИ 13-0100668, 14-01-00420.
функциями уравнений в частных производных, требуется сохранение непрерывности не только решения, но и его первой производной [2, 3,13, 17]. Требуемое свойство непрерывности производных искомого решения в аппроксимирующем численном решении на конечно-элементной сетке можно достаточно эффективно реализовать, используя интерполяционные эрмитовы элементы [17-19].
Данная мотивировка определила цель настоящей работы: построение вычислительных схем, алгоритмов и программ решения краевых задач методом конечных элементов, используя интерполяционные полиномы Эрмита (ИПЭ), сохраняющих непрерывность первых производных приближённого решения на конеч-ноэлементной сетке, реализация разработанных алгоритмов и программ в среде Maple-Fortran и анализ требуемых свойств аппроксимирующих решений в тестовых расчётах.
В работе представлены конечноэлементные вычислительные схемы решения краевых задач на примере задачи Штурма-Лиувилля с непрерывными, кусочно-непрерывными коэффициентными функциями, использующие ИПЭ. Выведены рекуррентные соотношения для вычисления в аналитическом виде ИПЭ с узлами произвольной кратности. Матрицы жёсткости и масс строятся в виде сумм интегралов от заданных коэффициентных функций исходного самосопряжённого дифференциального уравнения, вычисленных ИПЭ и их производных. Интегралы вычисляются с помощью гауссовых квадратур, а для кусочно-непрерывных полиномиальных коэффициентных и потенциальных функций — в аналитическом виде. Для аппроксимации таблично-заданных коэффициентных функций применимы интерполяционные формулы.
Построенные матрицы жёсткости и масс применяются для формулировки обобщённой алгебраической задачи на собственные значения, которая при малых размерностях матрицы решается, используя систему компьютерной алгебры. Эффективность - порядок аппроксимации по шагу конечноэлементной сетки и экономичность - время выполнения алгоритмов и программ в системе Maple для ленточных матриц размерностью порядка 100 и на языке Fortran при > 100 показана тестовыми расчётами точно решаемых задачах с непрерывными и кусочно-непрерывными потенциальными функциями.
Структура работы следующая. В разделе 2 дана формулировка задачи Штур-ма-Лиувилля и вариационного функционала. В разделе 3 получены рекуррентные соотношения для вычисления в аналитическом виде ИПЭ с узлами произвольной кратности. Представлена конечноэлементная схема с ИПЭ и описан алгоритм сведения краевых задач к алгебраическим. В разделе 4 представлен анализ эффективности и скорости сходимости вычислительных схем на тестовых расчётах. В Заключении обсуждаются результаты и перспективы применения построенных вычислительных схем.
2. Задача Штурма-Лиувилля и вариационный
функционал
Имеется самосопряжённое дифференциальное уравнение второго порядка относительно неизвестного решения Ф(^) на интервале х € Пх = (¿тт,%тах) [16]
1 д д
- 2Е) ФМ = о, о = —ш^АМш + Vм. (1)
Предполагается, что Д(х) > 0, %) > 0 — непрерывные коэффициентные функции, а V(х) — непрерывная потенциальная функция, которые имеют производные до порядка ктах — 1 > 1 в области г € Пг = [хшш, хтах], другие предположения для кусочно-непрерывных функций будут также рассмотрены в допустимых специальных случаях.
Для задачи Штурма-Лиувилля набор собственных функций Ф(-г) = Фт(%) € Н в пространстве Н, соответствующий набору собственных значений энергии
< ... ,< Ет < ..., удовлетворяет граничным условиям первого (I), второго (II) или третьего (III) рода при заданном значении параметра
(I) : Фт(г*) = 0, Ь = тт апё/ог тах, (2)
(II) : /1(2) ^ = 0, £ = тт апё/ог тах, (3)
(III) : ^^
с1г
Щгг)Фт(гг), £ = тт апё/ог тах (4)
и условию нормировки и ортогональности
г
(Фго(г)|Фго'(£)) = I ¡1 (г)(Фт(г))*Фт'(г)ёг = 5тт'. (5)
^тт
Решение краевых задач сводится к нахождению стационарных точек или минимальных значений вариационных функционалов [12,15]
Е(Ф, Е,гт[п,гтах) = у Ф*(г)(Б - 2Е)Ф(г)(1г = П(Ф, Е, гт[п, гтах)-
^тт
- /2(гтах)Ф*(гтах)Щгтах)Ф(гтах) + /2(гт[п)Ф*(гт[п)Щгт[п)Ф(гт[п), (6) где П(Ф, Е, , гшах) симметричный функционал:
П(Ф, Е,гт'ш,гтах)
М*) ™ ^ + ШФ*(г)У №)-
- Ш2ЕФ*(г)Ф(г)
йг. (7)
Здесь = 0 и ^ то для задачи дискретного спектра граничными усло-
виями второго рода и первого рода соответственно.
3. Вычислительные схемы МКЭ с ИПЭ
Вычислительная схема решения краевых задач (1)—(4) выводится из вариационного функционала (6), (7) на основе метода конечных элементов. В методе конечных элементов интервал [гтш,хтах] разбивается на подынтервалы, которые называются элементами. Размер этих элементов определяется достаточно гибко принимая во внимание физические свойства, качественное поведение и свойства гладкости искомого решения и его производных. Интервал А = [2:т1П^тах] покрывается набором п элементов А^ = хтах = 2™"], так что А = У™=1 А^. В результате имеем сетку
ф {г) [ХШ1П, ^тах] = {¿тт = гшш, гтах = + ^ = !,..., п - I,
*Тх = -С" + К = гтах], (8)
где = г^—х, ] = 2,... ,п точки сетки, величина шага = хтах - задаётся длиной элемента Ап.
г
3.1. Интерполяционные полиномы Эрмита
На каждом элементе Д,- задаём равномерную подсетку
Qhj(^max]
[ Zj
zr-\ = {Z(j-\)p = z^n, zU-i)p+r, r = l,...,p — 1, zjp = ^max}
с узловыми точками, zr = z^-i)p+r, определёнными по формуле ^■-Dp+r = ((р — г) zmin+г zmax)/P, r=o,...,p.
(9)
В качестве набора локальных функций {Nl(z, z™n, Zjmax)}\=0 , I
12 кг
на
r=0
каждом элементе г € [гтш, гтах] используем ИПЭ {{р^.(х)}Рг=0}к=0 1, заданные в узлах , г = 0,... ,р сетки (9). Значения функций р^(£) и их производных до порядка (Сах — 1), т.е. к = 0,..., Сах — 1, где ктах кратность узла гг, определяются по формулам [18]
к/ л Г Г dK <£r;(z)
(fir (^г') = Orr' Ок0, -
— farr' ^к
(10)
max
z=z
г
Для вычисления ИПЭ введём вспомогательную весовую функцию
-п max
p / \ К '
Wr (z) = П \JT~t) Г , Wr (Zr) = 1. (11)
Производные весовой функции представимы в виде произведения
= Wr (z) g К(z),
dKwr (z) к
dz к
где множитель вычисляется с помощью рекуррентных соотношений
ЯК*) = ^^^ + д1 (*) ^(г) (12)
с начальными условиями
^ 1 1 dwr (z) к
9o(z) = l, gr (z) = -г-)—— =
wr (z) dz , ' , z — zr'
r' = 0,r' = r
ИПЭ ^K(z) ищем в следующем виде:
ктах-1
Vt(z) = Wr (z) Y, агК> (Z — Zr )К . (13)
r
к =0
Продифференцировав (13) по z в точке zr и воспользовавшись (11), получаем
dK<pK(z)
d к
к
^-л к !
^ к''!(к' —
к"=0 у
с")!
дК -к
(z>r ^)ак
к''!.
(14)
2 = 2
Из (14) для коэффициентов а?'? получаем следующие выражения
_ 1 М? р? (х-)
а„
к! \ йг
'
? — 1 /1
- £ к/мV?(* )а?'К"к"! I . (15)
к"=0 Л )!
С учётом (10) из (15), получаем искомые коэффициенты разложения а. ИПЭ (13)
' 0, к' < к,
1/к'!, К = к,
к,к г
пк'к = < аг = \ к — 1
^ 1 /_ /' /' ^
- ¿^ (к — к')! Уг ' (гг)аК'К , к > к.
Р
Отметим, что степень всех ИПЭ не зависит от к и равна р = ^ ктах - 1.
г'=0
Далее будем рассматривать только ИПЭ с узлами одинаковой кратности кттах = ктах, г = 0,..., р. В этом случае степень полиномов равна р = ктах(р +1) - 1. Для таких полиномов введём следующее обозначение:
Нтхг+к(х, хтЫ, гтах) = <р'(г), г = 0,...,р, к = 0,...,ктах - 1. (16)
Такие ИПЭ образуют базис в пространстве полиномов степени р = ктах(р+1) -1 на элементе 2 € [гтт, гтах] с непрерывными производными до порядка ктах - 1 в граничных точках элемента.
При ктах = 1 ИПЭ совпадает с интерполяционным полиномом Лагранжа, который не имеет производных в граничных точках интервала г € [гтт, гтах]. ИПЭ при ктах = 1, 2,3 и р = 4 приведены на рис. 1.
2 = 2
к = к
Рис. 1. ИПЭ, совпадающий при кшах = 1 с интерполяционным полиномом Лагранжа (а) и ИПЭ при кшах = 2 (б, в) и кшах = 3 (г, д, е). Здесь р +1 = 5 -число узлов на элементе Д^ = [гШ1П = — 1,г,шах = 1]. Узлы сетки гг отмечены
вертикальными линиями
Видно, что значения ИПЭ ^та*р+к(г, хш1п, хгшах) и , х^", хЗ+х) (при г = р и г = 0) их производных до порядка кшах — 1 в точке хЗшах = г™" совпадают. Кроме того, граничные точки являются нулями кратностью кшах остальных ИПЭ независимо от длины элементов [хш1п, гшах] и [г™, ¿З+х]. Это обстоятельство позволяет построить базис кусочно-полиномиальных функций с непрерывными производными до порядка кшах — 1 на любом наборе А = Ц™=1 Аз = 1гшш, хшах] элементов Аз = [хш1п, хшах = гш+1].
3.2. Сведение задачи Шлурма—Лиувилля к алгебраической задаче
Построим дискретное представление решения Ф(-г) задачи (1)-(5), редуцированной к вариационному функционалу (6), (7) на конечноэлементной сетке,
П13 (г) [ *ш1п, = [ * 0 = ¿ш1п, = 1,...,пр — 1, гпр = гшах] (17)
с узлами XI = хзр = хушах = на сетке (г)[хшш, хшах], определённой в (8), и узловыми точками XI = Х(з-1)р+г, г = 0,...,р, подсеток О1 (г\хш1п, хЗшах], = 1,...,п, из уравнения (9). Решение Фн(х) ~ Ф(-г) ищем в виде разложения по базисным функциям N9(х) на интервале г € А = [хш1п, хшах]:
ФН(*) = Е Ф»N9(*), Фк(*) = Ф1т**, —( )
¡=0
Ф'lк,mzx + к,, (18)
г=гг
где Ь = (рп + 1)кшах — число базисных функций N¡9(х) и искомых коэффициентов
чения
С, ( г)[
Ф^ , которые при ^ = 1кшах + к есть значения производных к-го порядка функции Фн(х) в каждом узле 2 = Хк сетки ОЦ ( г)[хш1п,хшах], включая значение самой функции Фн(х) при к = 0.
Базисные функции N9 (х) — это кусочно-непрерывные полиномы порядка р' на интервале г € А = [хш1п, хшах]. Значения функции N¡9 (х) = N9^**+^) и её производных до порядка кшах — 1 равны нулю во всех узлах сетки (г), за исключением значения производной к-го порядка в узле XI равного единице, т.е.
^^У к.та* + к' (г)
йгк
= 5ц> б««', 1 = 0,...,пр, к = 0,...,кшах — 1.
Х = Х1
Для узлов х\ сетки (17), которые не совпадают с точками хшах сетки (8), т.е. при I = ]р, ] = 1,...,п — 1, полином N8(х) при ^ = (р(] — 1) + г)кшах + к определяется на элементе х € Аз как ИПЭ:
мд ( ) = [ Nкmaxr+к(г, ^ш1п, г,шах), х€ Аз; ( )
^(р(3-1)+г)кта*+к(г) = 0, х€ Аз. (19)
Поскольку точки хш1п и хшах есть нули кратности кшах, то такие кусочно-полиномиальные функции и их производные до порядка кшах — 1 являются непрерывными на всём интервале А. Соответствующие локальные функции — ИПЭ на рис. 1 показаны точечными, коротко-штриховыми и штрих-точечными линиями.
Для узлов XI сетки (17), совпадающих с одной из точек хшах сетки (8), принадлежащей двум элементам А и А +1 , = 1, . . . , п — 1, т.е. при = , кусочно-полиномиальная функция N!дкmaxj+K(x), производная которой порядка к в узле
= %тах равна единице, имеет вид:
= \ г, гтЦ, 4+Т), г £ А+1; (20)
0, х£ А ^ и А j+1.
Другими словами, кусочно-полиномиальные функции строятся как сшивка двух полиномов, первый из которых Мктахр+К(г, гтт, гтах) определён в области Аj, и значение его производной порядка к в узле ^ = хтах равно единице, а второй
Nк(z , zjja+l, - определён в области А j+l, и значение его производной порядка к в узле XI = хтах также равно единице. Эти кусочно-полиномиальные функции со всеми их производными до порядка ктах — 1 также непрерывны на интервале х £ А. Соответствующие локальные функции — ИПЭ на рис. 1 показаны сплошными и длинно-штриховыми линиями.
Подстановка разложения (18) в вариационный функционал (6), (7) сводит краевую задачу (1)-(5) к обобщённой алгебраической задаче относительно набора собственных значений Е собственных векторов Ф = {Ф^}1^—! :
(А — 2 Ен Б)Ф^ = 0. (21)
Здесь А = А + V + МШ1П — Мтах и Б симметричные матрицы жёсткости и масс размерностью Ь х Ь, где Ь = ктах(пр + 1):
Ак+1^+1 = £ 4; 12, (22)
0, к, I 2 )ео
= £ ^ 2, (23) и, к, I 2 )ео
В,1+1;,2+1 = £ , (24)
0, к, к)ео
^тах
. V dNll (г, гт1п, гтах) ^ (*, *т1п, *тах) ^ ^
Я;1 2 = / Ь &-^--^-Ь, (25)
~тт
К,12 = (г, гт1п, гтах)У(г, гт1п, гтах), (26)
~тт
Вк;12 = ЬШь (г, гт1п, ^^ (г, гт1п, гтах^г, (27)
~тт
где В = Ц £ {1,...,п}, Ь £ {0,...,р'}, 12 £ {0,...,Р= рктах(з — 1) + к, №2 = рктах(] — 1) +12}. Матрицы Мтах и Мтт размерностью Ь х Ь имеют только по одному ненулевому элементу: Мт1П = /2(гтш)Е(zmш) и Мт+х_ктах,ь+1_ктах = /2(гтах)К(гтах), соответственно.
1. Если коэффициентные /1 (г), /2(г) и/или потенциальная V(г) функции уравнения (1) заданы в табличной форме, то, используя подходящую интерполяцию
по ИПЭ, применяем интерполяционную квадратурную формулу, например, для матричных элементов Vii;i2 (zmm, zmax):
p ^max 1
Vli;l 2 = £ £ W^ 2;Kmaxr+KV ^ ( Z(j-1)p+r ),
r=0 к=0
где весовая функция Wivi2-i3 (zmin, zmax) даётся интегралом от произведения трёх ИПЭ:
max
Wit.i2.l3 = j h(z)Nh (z, zfn, zfax)Nh(z, zfn, zmax)Nh(z, zfn, zmax)dz.
Полученное выражение является точным для потенциальной функции V (г) в виде полинома степени меньшей, чем р. Иначе это разложение обеспечивает точность вычисления собственных функций и собственных значений до порядка р + 1.
Если интегралы (26) не вычисляются аналитически, то применяются квадратуры Гаусса [14,15] с р + 1 узлами на элементе, обеспечивающие теоретические оценки точности (29), достижимые в МКЭ,
Ai ^ ff ,dNh(z,zmm,zmax)
AL i2 =2^w9 h( Zg ) - • •
d
g=0
dNi2 (z, zfn, zmaax)
d
= Z9
VL 2 = £ Wg h( Zg )NLl ( Zg , Z?*, Zmax)V (Zg )NL2 (Zg , zf* , Z^) , (28)
g=0
K,2 = Y,wgh(zg)Nh (Zg, ^ *Tax)Ni2 (Zg, z^in, z^ax),
где гд = (р — д)гтт + дгтах и юд, д = 0,р — гауссовы узлы и веса ортогонального полинома порядка р + 1, определённого на элементе г € (гтш, хтах).
2. Расчёты в МКЭ обычно выполняются, используя переход к локальной координате г] € [— 1,1], связанной с абсолютной координатой г формулой: г = гтт+^(1+г1)/2, щ = Н^/2, где ^ = (гтах—гтт) — шаг конечноэлементной сетки (8). При этом локальные базисные функции и их производные пересчитываются по формулам
®(Z) = £ £ ^K.maXr+KiNKimaXr+Ki(r], -l, lW — )
г=0 к=0 4 '7
сСФ(х) у-л ^ йМктахг+к(г], —1,1) /¿г^ к
¿-^1 ¿-^1 к г+к Сп V Сг]
г=0 к=0 ' у '
3. Матрицы А^1ь2 , Вь1ь2 и VLlL2 - симметричные, размерностью ЬхЬ, где Ь = ктах(пр +1). Они состоят из п подматриц размерностью ктах(р+ 1) х ктах(р +1). Пересечения этих подматриц - блоки размерностью ктах х ктах. В матрицах и Аьгь2 внутри этих блоков есть нулевые элементы, их наличие связано с симметриями ИПЭ. В матрице VL1L■2 в общем случае внутри этих блоков нулевых элементов нет. Количество элементов во всех этих блоках каждой из матриц равно
z:-
(п(р2 + 2р) + 1)( ктах)2, и ширина ленты равна 2(рр + 1) —ктах. Примеры ленточной структуры матриц приведены на рис. 2.
Рис. 2. Структуры матриц Вь1ь2 и Аь1ь2 для шести (п = 6) элементов на всём интервале (хтш, гтах) в зависимости от кратности узлов ктах и числа подынтервалов р для схем седьмого порядка точности р' = 7. Слева направо:
(ктах,р) = (1, 7), (ктах,р) = (2, 3), (ктах,р) = (4,1), матрицы размерностью Ь х Ь, Ь = ктах(пр + 1), соответственно, равной 43 х 43, 38 х 38, 28 х 28, с общим числом элементов внутри блоков (п(р2 + 2р) + 1)(ктах)2 = 379, 364, 304 и шириной ленты 2(р' + 1) - ктах = 15,14,12
4. Матричные элементы (26) при учёте граничных условий второго рода (II) не изменяются. Для учёта граничных условий первого рода (I) в граничной точке zmiD из матрицы A вычёркивают первые строку и столбец, а для учёта граничных условий первого рода (I) в граничной точке zmax из матрицы A вычёркивают строку и столбец с номером L +1 — Kmax. Для учёта граничных условий третьего рода (III) в граничной точке zmin матричный элемент Ац даётся формулой: Ап = Ап + f2(zmiD)R(zmiD), а для учёта граничных условий третьего рода в граничной точке zmax матричный элемент Ai+i_Kmax i+1_Kmax даётся формулой AL+1_Kmax,L+1_Kmax = AL+1_Kmax,L + 1_Kmax — f2( ZmaX)R( ZmaX).
5. Решение обобщённой алгебраической задачи на собственные значения (21) для небольших размерностей до ~ 100 выполняется с помощью встроенных в пакет LinearAlgebra системы Maple процедур Eigenvectors и Eigenvalues. Для больших размерностей ~ 100 + 1000000 для решения задачи (21) используется метод итераций в подпространстве, реализованный на языке Fortran программой SSPACE [14], который эффективен для задач на собственные значения с симметричными ленточными матрицами большой размерности.
Теоретические оценки разности точного Фт(z) £ *Н2 и численного z) £ H к решений по норме H0 дают сходимость собственных значений и собственных функций порядка 2р' и р' + 1, соответственно [12]:
IEhm — Eml < C1 h2p', \\^>hm(z) — Фт(*)\\0 < C2hP' + 1, (29)
где h = max1<j<n hj наибольший шаг сетки.
4. Тестовые расчёты
Модифицированный потенциал Пёшля-Теллера. Рассмотрим решение задачи Штурма-Лиувилля для дифференциального уравнения Шрёдингера, в единицах Н = т = 1:
(Г+ WW — =
(30)
на всей оси z £ (—<х>, +то) с модифицированным потенциалом Пёшля-Теллера:
( ) 2 (cosh (a z))2, ( )
с параметрами a > 0 и Л > 0, для которого собственные функции и собственные значения известны в аналитическом виде. Это уравнение при V(z) = — k0n2(z) и 2Е = к2 = —332 применяется также для описания распространение электромагнитных волн в планарных тонкоплёночных оптических волноводах, которое устанавливает соответствие между профилем показателя преломления n( ) вол-новодного слоя и постоянной распространения 3. Здесь к0 = 2ж/Л0 - волновое число, где Ло- длина электромагнитной волны в вакууме [2]. Показатель преломления с зависимостью (31) от поперечной координаты применяется для расчёта мод в градиентных планарных волноводах [1].
Значения параметров Л = ll/2, a = l были выбраны таким образом, чтобы уравнение (30) с потенциалом V(z) = —99/8cosh(z)2 имело пять собственных значений 2 Ет = [—20, 25; —12, 25; —6,25; —2, 25; —0, 25].
Численные эксперименты на конечноэлементных сетках Q^ (z) [zmin = —40,
zmax = 40] c краевыми условиями (II), показали строгое соответствие теоретическим оценкам (29) для собственных значений и собственных функций. Для этого вычислялись коэффициенты Рунге
3 = l°g2
h h/2
h/2 h/4
1 = 1,2 (32)
на трёх сгущающихся сектах с абсолютными погрешностями собственных значений и собственных функций:
ah = \E^act — Et\, 4 = max,\ф™ас1 (z) — *™(z)\. (зз)
Из уравнений (33) получены теоретические оценки коэффициента Рунге сходимости собственных значений и собственных функций соответственно равные 2р' и (р' + 1).
В табл. 1 показаны коэффициенты Рунге для собственных значений и собственных функций третьего состояния, вычисленных с помощью схем до восьмого порядка р' = Kmax(p + 1) — 1 < 8 с различными Kmax и р, с шагом h = 0,125 для схем седьмого и восьмого порядка и h = 0, 0625 для остальных схем. В табл. 1 теоретические оценки коэффициента Рунге сходимости собственных значений и собственных функций соответственно равны 2р' и (р' + 1). Время расчёта Т(h = -1) пяти состояний дискретного спектра с шагом h = 1/32, число строк матриц Lh и число ненулевых элементов Nh представлено в трёх последних колонках.
Как видно из таблицы, для выбранного р' = 1 ^ 8 численные оценки коэффициента Рунге собственных значений отличаются от их теоретических оценок не более чем на 0,06 для р' = 1,..., 6 и не более чем на 0, 56 для р' = 7,8; для собственных функций — не более чем на 0, 2, т.е. строго соответствуют теоретическим оценкам (29). На рис. 3 показана зависимость абсолютных погрешностей ah = \e1xact — £h\ собственных значений ah = \e1xact — £h\ и собственных функций ah = maxze^h(zZ) \^xact(z) — xh(z)\ основного состояния в зависимости от шага сетки h конечноэлементных схем с ИПЭ для различных кmax и р. В двойной логарифмической шкале графики погрешности близки к прямым линиям с различным наклоном, что соответствует порядку аппроксимации р' = Kmax(p+ 1) — 1 при использовании ИПЭ с различными Kmax и р.
Таблица 1
Время расчёта Ти, число строк матриц Ьи и число ненулевых элементов при к = 1/32, коэффициенты Рунге для собственных значений и собственных функций третьего состояния (Е3 = -6, 25), вычисленных с помощью схем не более восьмого порядка р' = ктах(р +1) — 1 ^ 8 с различными ктах и р
_max fa p p Ru(E) 2p Ru^) p' + 1 Th Lh Nh
1 1 1 1,99 2 2,00 2 9,36 2561 7681
1 2 2 3,99 4 3,02 3 19,5 5121 20481
1 3 3 5,99 6 3,97 4 33,4 7681 38401
2 1 3 5,96 6 3,94 4 21,8 5122 30724
1 4 4 8,00 8 5,00 5 48,6 10241 61441
1 5 5 9,99 10 5,97 6 65,6 12801 89601
2 2 5 9,97 10 5,95 6 47,6 10242 81924
3 1 5 10,06 10 6,02 6 38,0 7683 69129
1 6 6 12,00 12 6,99 7 88,9 15361 122881
1 7 7 13,98 14 7,85 8 111, 17921 161281
2 3 7 13,87 14 7,77 8 82,3 15362 153604
4 1 7 13,57 14 7,59 8 59,6 10244 122896
1 8 8 15,99 16 9,09 9 139, 20481 204801
3 2 8 15,74 16 8,86 9 99,1 15363 184329
Рис. 3. Абсолютные погрешности of = IEfxact - Ev{ | и
|ФTact(z) - Ф\(г)1 собственного значения и собственной функции основного состояния задачи (30), (3) как функции от шага сетки h, вычисленные с ИПЭ при различной кратности узлов «max и числа подынтервалов р на элементе сетки
Для вычисления применялась программа KANTBP 1.1, реализованная на in Intel Fortran 77 с использованием типа данных QUADRUPLE PRECISION с 32 десятичными знаками. Время расчёта Т(h = 32) пяти состояний дискретного спектра с шагом h = 1/32 на компьютере 2 x Xeon 3,2 GHz, 4 GB RAM представлено в табл. 1. Схемы с фиксированным порядком р' имеют сходимость по шагу одного порядка, при этом время выполнения программы для схемы с большим значением ^max меньше, поскольку вычисления проводятся с матрицей меньшей размерности Lh x Lh с числом ненулевых элементов Nh, приведёнными в табл. 1.
Потенциал прямоугольной ямы. Для решения задачи в МКЭ с помощью ИПЭ в случае кусочно-непрерывных потенциалов в уравнении (30) применялась следующая техника аппроксимации (18) искомого решения, у которого первая производная непрерывна, а вторая производная терпит разрыв. Для кусочно-непрерывного потенциала V(z) = {Vi(z),z G (Cmin, Cmax)}, Cmax = Ci+f, где Vi(z)
являются р' + 1 кратно дифференцируемыми функциями, интервал (8), на котором определена задача, разбивался на подынтервалы [,гт1П, ,гтах] (,гтах = г™°) так, чтобы каждая точка 2 = С™ах = С™?, в которой решение задачи имеет разрывную вторую производную, совпадала с граничной точкой г^Х = ,гт1П двух соседних подынтервалов.
Рис. 4. Решения и их первые и вторые производные для основного состояния (сплошные кривые) и первого возбуждённого состояния (штриховые кривые) задачи с потенциалом прямоугольной потенциальной
ямы
Рассмотрим точно решаемую задачу Штурма-Лиувилля для уравнения (30) с потенциалом прямоугольной ямы шириной 2 а и высотой Vo/2, т.е.
2 V( z) = {Vo, N < а; 0,otherwise}. (34)
Чётные феп и нечётные собственные функции задачи (30) для такого потенциала имеют вид
С = К{е ^; * < -1; е"^ ^f^ , * € [-1,1]; е"^; * > 1},
cos(V-E - Vo)
С = N°{-e ^; * < -1; е"^ ^f—^l, * € [-1,1]; е"^; * > 1},
sin(VE - Vo)
где N® и N° - нормировочный множитель, который определяется из условия (5). При а =1, 2Vo = -50 задача Штурма-Лиувилля имеет пять собственных функций дискретного спектра (см. рис. 4)
2Ет = [-48,109146; -42,474904; -33, 232792; -20, 714111; -5, 965365].
Поскольку первые две собственные функции быстро убывают, то достаточно выполнить вычисления на интервале (г)[ zmin = -5, zmax = 5] c краевыми условием (II). Погрешность вычисления первых двух собственных значений представлена в табл. 2, а погрешность двух собственных функций и их первой производной на рис. 5.
Как видно из табл. 2, схемы с кратностями узлов ^max = 1 и ^max = 2 одинакового порядка точности р' = 3 и р' = 5 дают примерно одинаковую погрешность (при п = 20 погрешности порядка 10"2 и 4 ■ 10"6), а схема с ^max = 3 даёт существенно худшую погрешность (при п = 20 погрешности порядка 10"2). Это связано с тем, что вторая производная решения имеет разрывы в точках z = ± а.
В табл. 2 показаны коэффициенты Рунге (Ru) сходимости двух первых собственных значений из (32), вычисленных по схемам третьего р' = Kmax(p+1) -1 = 3 и пятого р' = nmax(p +1) - 1 = 5 порядков с различными ^max и р при h = 1/4, п = 40. Видно, что для выбранных схем третьего ' = 3 и пятого ' = 5 порядков численные оценки коэффициента Рунге отличаются не более чем на 0, 5 от их теоретических оценок, тогда как схема с ^max = 3, р = 1 даёт коэффициент Рунге схемы = 3 в три раза меньше теоретической оценки. Фактически эта схема
с = 3 = р' + 1 соответствует схеме второго порядка р' = 2, поскольку функция с разрывной второй производной аппроксимируется функцией с непрерывной второй производной.
Таблица 2
Абсолютные погрешности собственного значения энергии основного о±(Е0) и первого возбуждённого ^(Ег) состояний для потенциала прямоугольной
ямы (34) при а = 1 и 2У0 = -50. Коэффициенты Рунге (И,и) сходимости собственных значений из (32) при к = 1/4, п = 40 и его теоретические оценки 2р', даны в последних двух колонках
(кшах,р) р' (1,3) 3 (2,1) 3 (1,5) 5 (2,2) 5 (3,1) 5
1,93е-02 5,70е-02 2,47е-04 4,01е-04 1,48е-02
аГ/2(£о) 1,39е-03 3,15е-03 1,67е-06 2,59е-06 2,66е-03
/1=1/4/,-, N ^ ' (Ео) 4,44е-05 1,00е-04 3,82е-09 6,12е-09 3,51е-04
/=1/8/,-, \ ^ ' (Ео) 8,83е-07 2,21е-06 5,26е-12 8,59е-12 4,40е-05
а1=1/16(£0о) 1,48е-08 4,14е-08 2,22е-12 2,20е-13 5,50е-06
а?=1(^1) 9,96е-02 2,92е-01 6,44е-04 9,40е-04 6,70е-02
1=1/2/,-, \ ° 1 ' (Е1) 4,38е-03 1,14е-02 3,75е-06 5,66е-06 1,07е-02
=1/4 ° 1 ' №) 1,25е-04 3,08е-04 7,93е-09 1,27е-08 1,39е-03
=1/8 ° 1 ' №) 2,40е-06 6,33е-06 1,04е-11 1,74е-11 1,74е-04
=1/16 3,96е-08 1,14е-07 2,63е-12 2,06е-13 2,17е-05
М^о) Ии (£1) 2р' 5,65 5,70 6 5,50 5,60 6 10,3 9,99 10 9,51 9,53 10 2,99 3,01 10
Рис. 5. Разности между численными и точными значениями
.О^р 0'Р = Ф0 — Мг) для основного состояния (сплошные кривые)
.О^р {р = — для первого возбуждённого состояния (штриховые
кривые) для функций (сверху) и их первые производные (снизу) прямоугольной потенциальной ямы при п =10 элементов на интервале (—5, 5) для различной кратности узлов ктах и количества подынтервалов р на элементе сетки. Слева направо: (ктах,р) = (1, 3), (ктах,р) = (2,1),
(ктах,р) = (3,1)
Для иллюстрации этого факта на рис. 5 приведены погрешности и вычисления собственных функций и их первых производных. Как видно из рисунка,
среди схем третьего порядка р' = 3, схема с кшах = 2, р = 1 даёт в десять раз лучшее приближение для собственных функций и их первых производных, по сравнению со схемой с кшах = 1, р = 3. Видно, что производная погрешности аппроксимации для схемы с кшах = 1, р = 3 имеет разрывы в граничных точках элементов, которых нет для схемы с кшах = 2, р = 1. Однако схема пятого порядка р' = 5 с кшах = 3, р = 1 даёт худший результат даже по сравнению со схемами третьего порядка. При этом максимальные погрешности решения наблюдаются в окрестности границ ямы г = ± а, т.е. в этих точках подтверждается невозможность аппроксимации функции с разрывной второй производной функциями с помощью функций, имеющих непрерывные вторые производные.
5. Заключение
В работе представлены вычислительные схемы решения задач Штурма-Лиу-вилля методом конечных элементов с ИПЭ. Предложенные схемы порядка р' = Kmax(p+1) -1 при подходящем выборе числа р+1 узловых точек элемента сетки и кратности узлов Kmax сохраняют в приближённых решениях свойства непрерывности производных искомых решений. Показана эффективность вычислительных схем, реализованных в виде программ для матриц размерности 100 х 100 в системе Maple и большей размерности на языке Фортран в тестовых расчётах для точно решаемых квантовомеханических моделей с непрерывными и кусочно-непрерывными потенциальными функциями. Численный анализ приближенных решений в тестовых расчётах с непрерывными потенциалами показал, что порядок разработанных конечноэлементных схем, р' = Kmax(p +1) - 1, строго соответствует теоретическим оценкам. Выбор схемы большего порядка р' позволяет получить приближённое решение с более высокой точностью при большем шаге h конечноэлементной сетки при условии гладкости производной '-того порядка искомого решения.
Схемы с фиксированным порядком р' имеют сходимость по шагу одного порядка, при этом время выполнения программы для схемы с большим значением Kmax меньше, поскольку вычисления проводятся с матрицей меньшей размерности и с меньшей. Однако схемы с более высоким Kmax более чувствительны к погрешностям округления. Если производная порядка к искомого решения имеет точки разрыва, что имеет место в случае коэффициентных функций уравнения с разрывами производных порядка к - 2, то применение схем со значением кратности Kmax ^ к, т.е. ИПЭ с непрерывными производными большего порядка, не оправдано.
Для минимизации погрешностей предложенных конечноэлементных схем при заданных числе р + 1 узловых точек элемента сетки и кратности узлов Kmax состоит в выборе неравномерных подсеток. Также необходима разработка соответствующих конечноэлементных схем для случая, если вычисление коэффициентных или потенциальных функций занимает длительное время, при этом сетки, на которых вычисляются коэффициенты ОДУ, и её решения совпадать не будут.
Дальнейшее развитие и применение предложенных конечноэлементных схем, разработанных алгоритмов и программ решения задачи Штурма-Лиувилля связанно с анализом моделей молекулярных, атомных и ядерных систем, а также квантово-размерных систем, таких как квантовые точки, проволоки, и ямы в объёмных полупроводниках, и гладко-нерегулярных волноводных структур, которые описываются краевыми задачами для уравнения в частных производных с кусочно-непрерывными коэффициентными функциями и потенциалами.
Литература
1. Котляр В. В., Ковалев А. А., Налимов А. Г. Градиентные элементы микрооптики для достижения сверхразрешения // Компьютерная оптика. — 2009. — Т. 33. — С. 369-378.
2. Резанур Рахман К. М., Севастьянов Л. А. Задача одномерного рассеяния на ступенчатом потенциале с несовпадающими асимптотиками // Вестник РУДН. Серия «Физика». — 1997. — № 5. — С. 35-38.
3. Севастьянов Л. А., Егоров А. А. Теоретический анализ волноводного распространения электромагнитных волн в диэлектрических плавно-нерегулярных интегральных структурах // Оптика и спектроскопия. — 2008. — Т. 105. — С. 632-640.
4. Севастьянов Л. А, Егоров А. А, Севастьянов А. Л. Метод адиабатических мод в задачах плавно-нерегулярных открытых волноведущих структур // Ядерная физика. — 2013. — Т. 76. — С. 252-267.
5. Gusev A. A. Algorithm for Computing Wave Functions, Reflection and Transmission Matrices of the Multichannel Scattering Problem in the Adiabatic Representation using the Finite Element Method // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2014. — № 2. — С. 93-114.
6. Adiabatic Description of Nonspherical Quantum Dot Models / O. Gusev, A. A. Chuluunbaatar, S. I. Vinitsky, K. G. Dvoyan et al. // Physics of Atomic Nuclei. — 2012. — Vol. 75. — Pp. 1210-1226.
7. A Symbolic-Numerical Algorithm for Solving the Eigenvalue Problem for a Hydrogen Atom in the Magnetic Field: Cylindrical Coordinates / O. Chuluunbaatar, A. Gusev, V. Gerdt et al. // Lecture Notes in Computer Science. — 2007. — Vol. 4770. — Pp. 118-133.
8. Symbolic-Numeric Algorithms for Computer Analysis of Spheroidal Quantum Dot Models / A. A. Gusev, O. Chuluunbaatar, V. P. Gerdt et al. // Lecture Notes in computer Science. — 2010. — Vol. 6244. — Pp. 106-122.
9. Symbolic-Numerical Algorithms to Solve the Quantum Tunneling Problem for a Coupled Pair of Ions / A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar et al. // Lecture Notes in Computer Science. — 2011. — Vol. 6885. — Pp. 175-191.
10. Symbolic-Numerical Algorithm for Generating Cluster Eigenfunctions: Quantum Tunneling of Clusters Through Repulsive Barriers / S. Vinitsky, A. Gusev, O. Chuluunbaatar et al. // Lecture Notes in computer Science. — 2013. — Vol. 8136. — Pp. 427-440.
11. Models of Quantum Tunneling of a Diatomic Molecule Affected by Laser Pulses Through Repulsive Barriers / S. Vinitsky, A. Gusev, O. Chuluunbaatar et al. // Proc. SPIE. — 2014. — Vol. 9031. — P. 90311.
12. Strang G., Fix G. J. An Analysis of the Finite Element Method. — New York: Prentice-Hall, Englewood Cliffs, 1973.
13. Becker E. B., Carey G. F., Tinsley Oden J. Finite Elements. An Introduction. — New Jersey: Englewood Cliffs, Prentice Hall, 1981. — Vol. I.
14. Bathe K. J. Finite Element Procedures in Engineering Analysis. — New York: Englewood Cliffs, Prentice Hall, 1982.
15. KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach / O. Chuluunbaatar, A. A. Gusev, A. G. Abrashkevich et al. // Computer Physics Communications. — 2007. — Vol. 177. — Pp. 649-675.
16. ODPEVP: A program for Computing Eigenvalues and Eigenfunctions and their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined Sturm-Liouville Problem / O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich // Computer Physics Communications. — 2009. — Vol. 180. — Pp. 1358-1375.
17. Ramdas Ram-Mohan L. Finite Element and Boundary Element Applications in Quantum Mechanics. — New York: Oxford Univ. Press, 2002.
18. Березин И. С., Жидков Н. П. Методы вычислений. — М.: Физматлит, 1962. — Т. 1.
19. Самарский А. А, Гулин А. В. Численные методы. — М.: Наука, 1989.
UDC 517.958:530.145.6
Calculation Schemes for Solving Sturm—Liouville Problem by Finite-Element Method with Interpolating Hermite
Polynomials
A. A. Gusev*, L. L. Hait
* Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow region, Russian Federation, 141980
^ Belgorod State National Research University 85, Pobedy str., Belgorod, Russian Federation, 308015
Calculation schemes for solving Sturm-Liouville problem with first-, second- and third-type boundary conditions by finite-element method holding a continuity of derivatives of a required solution in its approximated solution are constructed. Recurrence relations for the calculation in analytical form of the interpolating Hermite polynomials with nodes of arbitrary multiplicity are derived. Using the interpolating Hermite polynomials, the basis piecewise-polynomial functions on finite-element grid with nonuniform step, approximating desired solution of the original problem are constructed and used for reduction to a generalized algebraic eigenvalue problem with banded stiffness and mass matrices. The stiffness and mass matrices are formed by sums of integrals containing the given coefficient and potential functions of the original self-adjoint second-order differential equation and the calculated interpolating Hermite polynomials and their derivatives on the finite element grid. The integrals are calculated using Gauss quadratures and in special cases, including the piecewise continuous polynomial coefficient and potential functions in analytical form. The efficiency and rate of convergence of the proposed calculation schemes and elaborated algorithms and programs implemented in Maple and Fortran is proved by benchmark calculations of exactly solvable Sturm-Liouville problems with continuous and piecewise continuous potential functions.
Key words and phrases: Sturm-Liouville problem, calculation scheme, finite element method, interpolation Hermite polynomials.
References
1. V. V. Kotlyar, A. A. Kovalev, A. G. Nalimov, Gradient-Index Elements of Microoptics for Superresolution, Computer Optics 33 (2009) 369-378, in Russian.
2. K. M. Rezanur Rakhman, L. A. Sevast'yanov, The One Dimensional Scattering Problem on Step Potential with Noncoinciding Asymptotics, Bulletin of PFUR, serie "Physics" (5) (1997) 35-38, in Russian.
3. L. A. Sevast'yanov, A. A. Egorov, Theoretical Analysis of the Waveguide Propagation of Electromagnetic Waves in Dielectric Smoothly-Irregular Integrated Structures, Optics and Spectroscopy 105 (2008) 576-584, in Russian.
4. L. A. Sevastianov, A. A. Egorov, A. L. Sevastyanov, Method of Adiabatic Modes in Studying Problems of Smoothly Irregular Open Waveguide Structures, Physics of Atomic Nuclei 76 (2013) 224-239, in Russian.
5. A. A. Gusev, Algorithm for Computing Wave Functions, Reflection and Transmission Matrices of the Multichannel Scattering Problem in the Adiabatic Representation Using the Finite Element Method, Bulletin of PFUR, series "Mathematics. Information Sciences. Physics" (2) (2014) 93-114.
6. O. Gusev, A. A. Chuluunbaatar, S. I. Vinitsky, K. G. Dvoyan, E. M. Kazaryan, H. A. Sarkisyan, V. L. Derbov, A. S. Klombotskaya, V. V. Serov, Adiabatic Description of Nonspherical Quantum Dot Models, Physics of Atomic Nuclei 75 (2012) 1210-1226.
7. O. Chuluunbaatar, A. Gusev, V. Gerdt, M. Kaschiev, V. Rostovtsev, V. Samoylov, T. Tupikova, S. Vinitsky, A Symbolic-Numerical Algorithm for Solving the Eigenvalue Problem for a Hydrogen Atom in the Magnetic Field: Cylindrical Coordinates, Lecture Notes in Computer Science 4770 (2007) 118-133.
8. A. A. Gusev, O. Chuluunbaatar, V. P. Gerdt, V. A. Rostovtsev, S. I. Vinitsky, V. L. Derbov, V. V. Serov, Symbolic-Numeric Algorithms for Computer Analysis of Spheroidal Quantum Dot Models, Lecture Notes in Computer Science 6244 (2010) 106-122.
9. A. A. Gusev, S. I. Vinitsky, O. Chuluunbaatar, V. P. Gerdt, V. A. Rostovtsev, Symbolic-Numerical Algorithms to Solve the Quantum Tunneling Problem for a Coupled Pair of Ions, Lecture Notes in Computer Science 6885 (2011) 175-191.
10. S. Vinitsky, A. Gusev, O. Chuluunbaatar, V. Rostovtsev, L. Le Hai, V. Derbov, P. Krassovitskiy, Symbolic-Numerical Algorithm for Generating Cluster Eigenfunctions: Quantum Tunneling of Clusters Through Repulsive Barriers, Lecture Notes in Computer Science 8136 (2013) 427-440.
11. S. Vinitsky, A. Gusev, O. Chuluunbaatar, L. Le Hai, V. Derbov, P. Krassovitskiy, Models of Quantum Tunneling of a Diatomic Molecule Affected by Laser Pulses Through Repulsive Barriers, Proc. SPIE 9031 (2014) 90311.
12. G. Strang, G. J. Fix, An Analysis of the Finite Element Method, Prentice-Hall, Englewood Cliffs, New York, 1973.
13. E. B. Becker, G. F. Carey, J. Tinsley Oden, Finite Elements. An Introduction, Vol. I, Englewood Cliffs, Prentice Hall, New Jersey, 1981.
14. K. J. Bathe, Finite Element Procedures in Engineering Analysis, Englewood Cliffs, Prentice Hall, New York, 1982.
15. O. Chuluunbaatar, A. A. Gusev, A. G. Abrashkevich, A. Amaya-Tapia, M. S. Kaschiev, S. Y. Larsen, S. I. Vinitsky, KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the CoupledChannel Hyperspherical Adiabatic Approach, Computer Physics Communications 177 (2007) 649-675.
16. O. Chuluunbaatar, A. A. Gusev, S. I. Vinitsky, A. G. Abrashkevich, ODPEVP: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined Sturm-Liouville Problem, Computer Physics Communications 180 (2009) 1358-1375.
17. L. Ramdas Ram-Mohan, Finite Element and Boundary Element Applications in Quantum Mechanics, Oxford Univ. Press, New York, 2002.
18. I. S. Berezin, N. P. Zhidkov, Computational Methods, Vol. I, Pergamon Press, Oxford, 1965.
19. A. A. Samarski, A. V. Gulin, Numerical Methods, Nauka, Moscow, 1989, in Russian.
