CC BY
72
Математическое моделирование
УДК 519.633.2, 517.958, 530.145.6 Алгоритм численного решения параметрической
задачи Ш^турма-Лиувилля и вычисления производных от решения по параметру методом
конечных элементов
О. Чулуунбаатар
Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, Московской обл., г. Дубна, Россия, 141980
Представлен экономичный алгоритм численного решения параметрической задачи Штурма-Лиувилля с краевыми условиями третьего рода на конечном интервале методом конечных элементов с автоматическим выбором сдвига спектра. Алгоритм включает вычисление с заданной точностью набора ~ 10—50 собственных значений, собственных функций и их первых производных по параметру, а также интегралов — матричных элементов между собственными функциями и их производными по параметру. Эффективность алгоритма продемонстрирована численным анализом интегрируемой задачи с параметрическими краевыми условиями третьего рода.
Ключевые слова: параметрическая задача Штурма-Лиувилля, метод конечных элементов, метод Канторовича.
1. Введение
Математические модели физических процессов фотоионизации и лазерно-стимулированной рекомбинации атома водорода в однородном магнитном поле при воздействии лазерного излучения, осевого каналирования одноимённо-заряженных частиц в канале кристалла, возбуждения и де-возбуждения волнового пакета атома водорода в однородном магнитном поле при воздействии последовательности сверхкоротких лазерных импульсов изучались в недавних работах [1—4]. Подобные задачи возникают также при изучении модели процессов фотоабсорбции в квантовой яме с водородоподобным примесным центром [5].
Математические модели этих процессов объединены объектом численного исследования, которым являются сингулярные спектральные, краевые и эволюционные задачи для многомерного уравнения Шрёдингера в конфигурационном пространстве [6,7]. Наличие нескольких потенциалов взаимодействия между частицами или частиц со внешними полями приводит к разбиению конфигурационного пространства на подобласти, в каждой из которых доминирует тот или иной потенциал. Кроме того, сингулярная краевая задача с помощью асимптотических разложений решения в асимптотической области конфигурационного пространства редуцируется к регулярной переносом асимптотических краевых условий в виде условий третьего рода на границу конечной области [8].
Для решения краевой задачи в сложной области проекционными методами строятся базисные функции с нелинейными параметрами в одной из подобластей или в каждой подобласти. На границах этих подобластей для согласования составных базисных функций в общем случае следует использовать условия третьего
Статья поступила в редакцию 6 февраля 2009 г.
Работа выполнена в рамках темы ОИЯИ «Математическая поддержка экспериментальных и теоретических исследований, проводимых ОИЯИ 05-6-1060-2005/2010» и поддержана РФФИ (грант 08-01-00604-а «Математическое моделирование динамики лёгких атомов и молекул под действием быстрых частиц, лазерных импульсов и магнитных полей»).
Автор благодарит С.И. Виницкого, А.А. Гусева, И.В. Пузынина и Л.А. Севастьянова за сотрудничество и поддержку.
рода. При этом количество нелинейных параметров в вариационном функционале существенно увеличивается, что требует для их оптимизации значительных ресурсов ЭВМ.
Метод Канторовича (МК) [9] даёт возможность построить экономичный алгоритм вычисления параметрических базисных функций, которые учитывают все особенности сложных областей и условия согласования между ними. Использовании такого базиса позволяет понизить размерность исходной краевой задачи сведением к системе бесконечного числа обыкновенных дифференциальных уравнений с матрицами переменных коэффициентов и краевыми условиями третьего рода [6,7]. При этом требуется исследовать скорость сходимости разложения искомого решения по числу базисных функций. Дискретизация краевых задач методом конечных элементов (МКЭ) [10,11] приводит к последовательности параметрических алгебраических задач, для решения которых требуется разработка экономичных алгоритмов и программ [12,13]. Требования заданной точности результатов вычислений и экономичности алгоритмов определяются необходимостью вычисления набора элементов для матриц переменных коэффициентов размерностью порядка ~ 10-50. Заметим, что при вычислении изолированного собственного решения, зависящего от параметра, методами Ньютона или обратных итераций в качестве начального приближения или сдвига спектра обычно используют собственные значения и их производные, вычисленные при достаточно близком предыдущем значении параметра. В случае плохо изолированных решений такой подход требует проведения вычислений на достаточно плотной сетке значений параметра из заданного интервала с учётом дополнительных условий ортогональности вычисленного набора решений [14]. Это затруднение преодолевается в рамках метода обратных итераций в подпространстве с автоматическим выполнением условия ортогональности [11].
В данной работе представлен экономичный алгоритм численного решения с заданной точностью набора ~ 10-50 собственных значений, собственных функций и их первых производных по параметру параметрической задачи Штурма-Лиувилля с краевыми условиями третьего рода на конечном интервале. МКЭ применяется для дискретизации задачи и построения численных схем вычисления собственных функций и их производных по параметру с точностью порядка 0(кр+1), а также собственных значений, их производных и матричных элементов 0(к2р) по шагу И конечно-элементной сетки [10,11]. Выбор порядка аппроксимации р зависит от гладкости искомого решения. Представлен и апробирован алгоритм для вычисления нижней оценки наименьшего собственного значения обобщённой параметрической алгебраической задачи на собственные значения. Этот алгоритм позволяет автоматически определить сдвиг спектра и сэкономить число итераций при вычислении набора собственных решений методом обратных итераций в подпространстве с помощью подпрограммы ББРАСЕ [11].
2. Постановка задачи
Рассмотрим параметрическую задачу Штурма-Лиувилля на конечном интервале X £ Пг = (Zшln, %шах):
Ь(г; р)ф(г; р) = J1^dLf2(z)dL + П(р,г)^1 ^ р) = е(р)^(г] р)' (1)
Здесь р £ Пр = [рш[п, ршах] — вещественный параметр, е(р) — собственное значение, зависящее от параметра р. Предполагается, что функции /т_(г) > 0, > 0, и и(р,г) непрерывны на интервале г £ Пг, а функция и(р,г) дифференцируема по параметру р. Собственная функция ф(г; р), зависящая от параметра р, подчиняется краевым условиям третьего рода на левом и правом концах интервала г £ ° г, т.е. на границе д области = ° дПг:
11(р,р1,Х1(р),ф(г; р)) = Р) + Х1(р)ф(г; р) = 0 г = гшп (2)
12(р,р2,А2(р),ф(г;р)) = р2/2+ \2{р)'ф{г; р) = г = гтах, (3)
где Ах(р), А2(р) — дифференцируемые по параметру р функции. Далее считаем, что = 0 для краевого условия первого рода, и /лj = 1 для остальных случаев.1 В представленном ниже алгоритме сначала вычисляется набор ,]тах наименьших собственных значений ех(р) < е2(р) < ... < е^тах(р), причём ех(р) ^ а(р),
и соответствующих им собственных функций {ф^ (х;р) }^=а1х удовлетворяющих условиям ортогональности и нормировки
^■т ах
У р)ф(х;р)Ах = 8ц, (4)
^т 1п
где — символ Кронекера и а(р) > —<х — нижняя оценка наименьшего собственного значения 1( р). Далее вычисляется набор производных собственного значения де^ (р)/др и собственной функции дгф^ (х;р)/др.
Дифференцирование краевой задачи (1)-(3) по параметру р приводит к решению неоднородной краевой задачи относительно искомых параметрических производных собственной функции д'фj (г; р)/др и собственного значения де^ (р)/др
/л ч / ЛЛдфз (*;Р) (ди (Р, д€з (Р)\ I ( \ ( Дг;р) - ^ (р))-—-= - ^^р---о^) '(z;Р), (5)
0, (6)
dl i(p,p,\,Xi(p),^j (zm[n;p))
dp
■?/]-( ^ • л П
0. (7)
dl 2 (p,p2,\2(p),^j (z max) p)
д р
Умножая уравнение (5) слева на собственную функцию фъ(г;р) и интегрируя на интервале г е , получаем соотношения
^т ах
| Мг)фг(г;р)(Ь(г;р) - Ч(р)) dz =
= (£i(p) - 4 (p)) I fl(z)lpi( Z-p)-d-dp-+ fij (p, Zmn Zmax) =
Z-
min
^m ax
/1(;-,р)(- Цр>)ф,(8)
-^тт
где функция ¡^ (р, ^тт, ^тах ) в зависимости от выбранного типа краевых условий задана в следующем виде
Iij (р, Zmm, ^тах)
-'-Если на левом конце интервала /2(2) = (z — zmin)n~K(z), -K(zmin) > 0, п ^ 1, а функция fi(zmin + 0)U(p,zmin + 0) либо ограничена, либо стремится к то собственная функция подчиняется
краевому условию типа (2) при Ai(p) = 0
lim f2(z)--- = 0,
z^zmin + ° dZ
которое есть следствие естественной ограниченности искомого решения |^(zmin; р) | < те [8]. Аналогичное краевое условие типа (3) задаётся на правом конце интервала при соответствующих требованиях к функциям fi(z), /2(2) и U(p,z).
Z
др
Р.)' (гшах;р)
- ^дрг^Фг&шт, р)' (гшт',р), если р1Х1(р) = 0, Р2Х2 (р) = 0,
= < дЧрг'( (гша*;р), если р1Х1(р) = 0, (р) = 0, (9)
-д15рр)'г(р)'о (*шт;р), если р^1(р) = 0, Ц2Х2 (р) = 0, 0, если р1Х1(р) = 0, р2Х2 (р) = 0.
В силу этих соотношений получаем формулу для вычисления параметрической производной собственного значения
^шах
= / (г;р)-Щр^-' (г;р)йг + ¡п(р, , ^х). (10)
-^шт
Поскольку оператор Ь(г; р) — е, (р) вырожденный, то задача (5)-(7) имеет бесконечное множество решений. Из условия нормировки (4) следует требуемое дополнительное условие, позволяющее выделить единственное решение задачи (5)-
(7)
^шах
I ШФ, (г;р) <Ъ = 0. (11)
-^шт
Данная постановка задачи (5)—(11) связана с вычислением соответствующих интегралов — матричных элементов, используемых при решении многомерных краевых задач в рамках МК [6,7]:
Яг,(р) = — ¡1 ,
■^-шт ^-шах
д'з (г;р), - р
(12)
д'г(г;р) д',(г; р) - р - р
Нц(р)= к (г^Т^^Т^йг.
Чтобы избежать явного вычисления производных решений по параметру, для вычисления набора матричных элементов {(^гз(р)}1ш=1 и {—^(р)}^1^, на первый взгляд, можно использовать следующие формулы:
2-т ах
е Ар) — ^ (р) ^ др
шт
р, ¿шт, ¿шах) ■ / ■ гл ( \ п /ю\
+ ег(р)—Ч (р) , % =(гг (р) = 0, (13)
ах
Нгз(р) = — Шц (р). (14)
1 = 1
Однако разложение матричных элементов {Нг,(р)}\™=1 по формуле (14) имеет медленную скорость сходимости по числу матричных элементов ((ц(р), т.е., с увеличением размерности )х для достижения заданной точности требуется 1тах ^ ]ш)х [12,13]. Это обстоятельство приводит к существенным затратам ресурсов ЭВМ и определяет цель данной работы: разработать и апробировать экономичный алгоритм решения с заданной точностью параметрических краевых задач типа (1)—(11) и вычисления с той же точностью производных по параметру от решений, а также набора соответствующих интегралов (12).
0
Поскольку исходная задача (1)—(3) линейная, то для обеспечения непрерывности набора решений и матричных элементов по параметру на интервале р С необходимо использовать дополнительное условие на границе интервала интегрирования, например фj (г; р) > 0 в окрестности точки г = гшах.
3. Формулировка алгебраических задач МКЭ
Сформулируем численные алгоритмы для вычисления решений параметрической краевой задачи (1)—(3) и их производных по параметру р. Вычислительные схемы построены в рамках МКЭ на основе вариационного функционала Рэлея-Ритца
{•'ша^ 2 \
/ (^К ¿ь) + ш)и(р,*№2(*;р)У* +
-1
+ 9(р1 ^'ш\n, ^шах)
[ х I I (г; I ' (15)
где функция д(р, ^шт , ^шах ) в зависимости от выбранного типа краевых условий (9) задана в следующем виде:
9 (р? Zшm, ^шах)
\2(р)Ф(¿шах; р)Ф(¿шах; р)
-\1(р)ф(гтт;р)ф( ¿ш1п;р), если 1ц\1(р) = 0,
= < \2(р)Ф( ¿шах;рЖ ¿шах; р), если щ\1(р) = 0, -\1(р)ф(гшт;р)ф( ¿ш1п;р), если 1ц\1(р) = 0,
0,
если ц1Х1(р) = 0,
Р'2^2 (р) = 0, Ц2М(р) = 0, Ц2М(р) = 0, Ц2^2(р) = 0.
(16)
МКЭ решения краевых задач для дифференциального уравнения второго порядка состоит в том, что конечный интервал разбивается на подынтервалы ^з = [гз — 1, ] = и^=1 ^з), называемые конечными элементами [10]. В каждом из подынтервалов ^^ определяются узлы2 г? г и строятся интерполяционные полиномы Лагранжа
Ф1г(г)
г=0,г=г ^ 3,г /^3,г'
(17)
2 Мы используем равномерное разбиение каждого из подынтервалов
hj -
— 1 +--г, hj = г^ — ¿2-1, г = 0,р.
р
z
С помощью этих полиномов фР г (х) определяется набор локальных функ-
ций N (р):
*Г(*) =
ф10(г), ге А1,
0, хе Д1,
фрьг(х), хе Д,,
Ф1Р(*),
фj+l,o(z), 0,
ФPn,p(z), 0,
х е д, ,
X е д , х е д,+и
хе д,-и
х е Дп, х е дп,
= 0,
I = Г + р(] — 1), г = 1,р — 1,
(18)
l = jp, 3 = 1,п — 1
I = пр.
Кусочно-дифференцируемые функции {^(х)}^=0, Ь = пр, формируют неортогональный базис в пространстве полиномов порядка р. Функция ф(х; р) е Трр ~ х ^р) аппроксимируется конечной суммой локальных функций ^(х)
Ф(х;р) = ,р)Щ (х),
(19)
1=0
т.е. функция ф(х; р) имеет лишь обобщённые частные производные первого порядка и принадлежит пространству Соболева W1 [10]. В результате подстановки выражений (19) в вариационный функционал (15) и его минимизации [10,11] получаем обобщённую алгебраическую задачу на собственные значения для определения численных решений ер и -фр
Ар'фр = е рВр'фр, Ар = Ар + М.
(20)
Здесь М — диагональная матрица с нулевыми элементами, кроме первого и последнего элементов, которые равны Х2(р) (или нулю) и —Х1(р) (или нулю), соответственно. Матрицы жёсткости Ар и массы Вр — симметричные и ленточные матрицы (с шириной ленты 2р + 1 каждая), причём Вр — положительно определена. Эти матрицы представимы в виде
Ар = ^ ар, Вр = ^ Ь
(21)
3=1
3 = 1
Локальные матрицы ар и Ьр вычисляются по формулам
+1
(ар)яг =
1
+1
(Ь,р)
р\дг _
Ь(хЩА (х)фр,г (г)2 ^,
1
(22)
х = х,-1 + 0,5Из (1 + 77), д,г = 0,р.
Интегралы (22) вычисляются с помощью гауссовской квадратуры порядка р + 1
ю
* Г лАФ1ч (га )ЩГ (г9)
1 /2 {га
а=о ^
^ ^ „ + МХа)П{р^а{,аЛ|щ,
Ь2 ¿г
¿г
{ЦГ = Е ^а Ж, (*а Щ,Г (га) ^ ™а,
а=о
0
где Zg = Zj-\ + 0, 5hj(1 + r]g), rjg и wg, g = 0,p — гауссовские узлы и веса [10]. Это обстоятельство позволяет решать параметрическую краевую задачу (1)—(3) на интервале z Е Ùz с учётом естественной ограниченности искомого решения1.
Для численного решения имеют место следующие оценки погрешности [10] для положительно определённой матрицы Ар
|е,(p) - €j \ < ci|ej(p)\h2p, ||ф(Z]p) -ф"||о < c2|(p)\hp+1, (24)
zmax ,
где ||v(z; p)||0 = f f1(z)v2(z;p)dz; €j(p), фj(z\p) — точные решения; ej, ф" со-
zmin
ответствующие численные решения; h — максимальный шаг сетки; m — номер решения; ci и с2 — положительные константы, не зависящие от шага h.
3.1. Алгоритм вычисления параметрических производных от решения и матричных элементов
С помощью представленной выше процедуры дискретизации неоднородная краевая задача (5)—(7) сводится к системе линейных неоднородных алгебраических уравнений
L a-f = (AP - fhBP) = „, b = - ( - BP ) ф". (25)
dp dp Условие нормировки (4) и дополнительное условие (11) принимают вид
(фн)Т Bрф" = 1, Bрф" = 0. (26)
В результате параметрическая производная собственного значения (р), матричные элементы Н^(р) и Qij(р) вычисляются по формулам
de'j(p) = (ф"\г dAV о" (p) = (ФАТВр?ф" Hj (p)=(ж)ВР-
% (p)=[ Вр
Поскольку е н является собственным значением задачи (20), матрица Ь в (25) вырожденная (или близка к вырожденной в случае численного решения). Соответствующий алгоритм решения задачи (25), (26) реализуется в виде последовательности трёх шагов:
ТТТаг 1. Вычисление решений V и w двух вспомогательных систем линейных неоднородных алгебраических уравнений
IV = Ь, Lw = ^ (28)
с невырожденной матрицей Ь и правым частями Ь и d
г =( LSS', (s -S )(s> -S ) = 0, Lss' \ 5aa,, (s -S )(s ' -S ) = 0,
bs, S = S, л _ f Lss, S = S,
j
h =i bs, S = S, И = f
°s = \ 0, s = S, as =\ 0, s = S где S — номер максимального по модулю элемента вектора Врф
ТТТаг 2. Вычисление трёх вспомогательных коэффициентов 71, 72 и D^
7 =
71
(Ds - 72)'
71 = v
т
72 = w
T W^h, Ds = (Bp фh)s.
(30)
Т!аг 3. Вычисление вектора д "ф /др
J vs
7,
vs — 7WS, s = S,
др 1 7, s = S.
(31)
Отметим, что вторая оценка (24) также имеет место для решения (5)—(11). Этот факт гарантирует такую же точность для производных собственных значений и собственных функций, а также матричных элементов в рамках представленного алгоритма, т.е. имеют место следующие оценки для положительно определённой матрицы Ар
де 0(р) Оеh
д р д р
< сз| Ч Ш2р,
дф (Z]p) д^
д р
д р
< Iej (pW+1, (32)
Qij(р) — Qh0(Р) < 1ег(р)е0(p)Ih2p, Нго(р) — Щ(р) < се|€i(р)е0(р)^, (33) где Сз, С4, С5 и се — положительные константы не зависящие от шага h.
0
3.2. Алгоритм для вычисления нижней оценки наименьшего собственного значения
В общем случае краевая задача (1)-(3) имеет как отрицательные, так и положительные собственные значения, поскольку мы не требовали выполнения условий U(р, z) ^ 0, Х2(р) ^ 0 и \1(р) ^ 0, обеспечивающих не отрицательность собственных значений 0 ^ е1(р) < е2(р) < ... < ejmax (р). В силу этого обстоятельства матрица жёсткости Ap = Ap(р) в общем случае не является положительно определённой.
Поэтому затруднительно получить достаточно близкую нижнюю оценку а = а(р) < efh(р) наименьшего собственного значения е1(р) для задачи (20), например, используя известные теоретические оценки кругов Гершгорина [15],3 при каждом вещественном значении параметра р из некоторого заданного интервала Qp. На основе известного свойства разложения F = LDLT для симметричных матриц F: для положительно определённых симметричных матриц F все элементы диагональной матрицы D положительны [16], нами предложен и апробирован следующий алгоритм:
ТШаг 1. Выполнить LDLT факторизацию для Ap — aBp.
ТТТаг 2. Если некоторые элементы диагональной матрицы D меньше чем 0, тогда положить а = а — 1 и перейти к ТШагу 3, иначе перейти к ТТТагу 5.
ТШаг 3. Выполнить LDLT факторизацию для Ap — aBp.
ТТТаг 4. Если некоторые элементы диагональной матрицы D меньше чем 0, тогда положить а = а — 1 и перейти к Тагу 3, иначе положить а = а — 0.5 и перейти к ТТТагу 8.
ТТТаг 5. Положить а = а + 1 и выполнить факторизацию LDLT для Ap — аВ1'.
ТТТаг 6. Если все элементы диагональной матрицы D больше чем 0, тогда положить а = а + 1 и перейти к Тагу 5.
ТТаг 7. Положить а = а — 1.5.
ТТТаг 8. Вывод искомого значения а* = а.
Начальное приближение а(р) = а0 в алгоритме (1-8) задаётся a priori, например, ао = 0 или с помощью известных асимптотических оценок ао < ^л(р) [1].
3или а{р) = infze[zmin,zmax] и{p,z) при \2(р) > о и \1{р) < 0 для задачи (1)-(3).
Выполняя данный алгоритм (1—8), получаем нижнюю оценку наименьшего собственного значения задачи (20), и всегда е1(р) — а* ^ 1.5, т.е. полученная матрица Ар—а*Вр положительно определена. В результате (20) сводится к задаче
(Ар — а*Вр) ф} = ен Врф}, ен = ен + а*, (34)
которая имеет положительные собственные значения 0 < ¿^(р) < ... < (р).
Таким образом, алгоритм (1—8) позволяет автоматически определить сдвиг а* спектра задачи (20) и сэкономить число итераций при вычислении набора решений методом обратных итераций в подпространстве с помощью подпрограммы ББРАСЕ [11] при любом значении параметра р из заданного интервала Пр. Для решения последовательности краевых задач (1)—(3) на сетке П}ц [ршш,ршах] = {р1 = Pmin,Pi+1 = Pi + Ы, р^х = ртах} с достаточно мелким переменным шагом hi значения а(р^{) можно оценить по формуле а(р^{) = + (д£\(р^/др — \де}}(р{)/др\ — Щ)Ы/2, & > 1, вычисленные на предыдущем шаге рi, уже не используя алгоритм (1—8). В этом случае алгоритм (1—8) используется, например, только при значении р = рт;п или р = ртах.
Указанные выше алгоритмы были реализованы в виде комплекса программ РАБОБР [17] на языке Фортран 77.
4. Точнорешаемая параметрическая модель
В качестве тестового примера использовалась задача Штурма-Лиувилля
д(в;р)
?в2
= ej (p)iPj ( 6;p) (35)
с краевыми условиями по угловой переменной на дуге —ж/6 ^ в ^ ж/6 сектора круга радиуса р, который рассматривается здесь как параметр:
«M Tpc^j (0;p) = 0, 0 = Т6. (36)
Для краевой задачи (35)—(36) известны чётные и нечётные собственные функции и собственные значения через решения трансцендентного уравнения, которые вычисляются с компьютерной точностью, а также аналитические выражения для набора матричных элементов {Oij(p)}'jmj=1 и {Hij(p)}3i'm=X1, которые выражаются только через соответствующие собственные значения [18]. Ниже приведены результаты численного анализа только для чётных функций с краевыми условиями по угловой переменной -ж/6 ^ в ^ 0 при значениях параметров с = -1, к = ж/6, ^i(p) = -pcK и M(p) = 0:
ММ _pcmPj(S;p)=0, в= -6, ММ = 0, в = 0. (37)
4.1. Численные результаты и перспективы
При вычислении собственных значений была выбрана конечно-элементная сетка П}в [$тт, $тах] = {$тт = —ж/6(200) 0тах = 0} с элементами четвёртого порядка (р = 4). Число в скобках обозначает число конечных элементов на по-динтервале. Расчёты выполнялись при значениях радиуса р, заданных на сетке П}р[ртт,ртах], применяемой для решения соответствующих краевых задач для системы обыкновенных уравнений второго порядка с переменными коэффициентами по радиальной переменной р, полученных при редукции исходной двумерной краевой задачи МК в параметрическом базисе (35), (36) [7,18]. Наборы вычисленных собственных значений {(р), матричных элементов ^^(р)}^^ и
Рис. 1. Вычисление собственных значений е^ (р) и матричных элементов Qij(p) и Нг ¡(р) в зависимости от параметра р при с = —1 и jmax = 6
{Hij— переменных коэффициентов радиальных уравнений в зависимости от параметра р при с= — 1 и jmax = 6 приведены на рис. 1.
Численные эксперименты показали строгое соответствие теоретическим оценкам погрешности для собственных значений, собственных функций (24) и их производных по параметру (33). В частности, были вычислены значения коэффициентов Рунге [7] на трёх вдвое сгущающихся сетках для собственных значений, их производных и матричных элементов в пределах 7.99 ^ 8.01, а для собственных функций и их производных в пределах 4.99 ^ 5.01, что соответствует теоретическим оценкам погрешности выбранной аппроксимации порядка р = 4. В табл. 1-3 приведено сравнение с аналитическими результатами для собственных значений, их производных по параметру и матричных элементов в точке р = 2. Расчёты выполнялись с заданной точности 10-15 на компьютере 2 x Xeon 3.2 GHz, 4 GB RAM, используя Intel Fortran 77 и тип данных real*16, обеспечивающий 33 значащие цифры. Время счета для данного примера составляет 2 секунды.
Таблица 1
Относительные погрешности между точными и вычисленными собственными значениями ех = (р) — £^(р) 1/1€] (р)| и их производными
£2 = \дej (р)/др — де1;:(р)/др\/\е.} (р)| при р = 2 и = 6
3 £1 £2
1 0.9490E-27 0.2149E-26
2 0.1018E-21 0.2628E-22
3 0.3431E-19 0.1974E-20
4 0.9222E-18 0.2313E-19
5 0.9364E-17 0.1458E-18
6 0.5623E-16 0.8309E-15
Таблица 2
Абсолютные погрешности между точными и вычисленными матричными элементами (р) — Н^ (р)| при р = 2 и ]ш&х = 6
i /3 1 2 3 4 5 6
1 0.406E-26 0.113E-23 0.507E-22 0.526E-21 0.124E-19 0.307E-16
2 0.113E-23 0.191E-23 0.181E-22 0.689E-22 0.188E-20 0.269E-16
3 0.507E-22 0.181E-22 0.661E-22 0.187E-21 0.508E-20 0.118E-16
4 0.526E-21 0.689E-22 0.187E-21 0.634E-21 0.577E-20 0.922E-18
5 0.124E-19 0.188E-20 0.508E-20 0.577E-20 0.107E-19 0.924E-17
6 0.307E-16 0.269E-16 0.118E-16 0.922E-18 0.924E-17 0.392E-16
Таблица 3
Абсолютные погрешности между точными и вычисленными матричными элементами (р) - . (р)| при р = 2 и ]тах = 6
i/3 1 2 3 4 5 6
1 0.138E-33 0.550E-23 0.419E-21 0.492E-20 0.856E-21 0.646E-15
2 0.550E-23 0.359E-35 0.460E-21 0.578E-20 0.226E-18 0.528E-15
3 0.419E-21 0.460E-21 0.165E-35 0.336E-20 0.151E-18 0.803E-16
4 0.492E-20 0.578E-20 0.335E-20 0.355E-36 0.887E-19 0.132E-15
5 0.205E-19 0.242E-19 0.209E-19 0.512E-20 0.555E-36 0.151E-14
6 0.349E-16 0.417E-16 0.417E-16 0.424E-16 0.434E-16 0.712E-36
5. Заключение и обсуждение результатов
Представлен экономичный алгоритм численного решения параметрической задачи Штурма-Лиувилля с краевыми условиями третьего рода на конечном интервале МКЭ. Алгоритм включает вычисления с заданной точностью набора ~ 10-50 собственных значений, собственных функций и их первых производных по параметру и интегралов — матричных элементов между собственными функциями и их производными по параметру. Показано, что дискретизация задачи МКЭ, разработанные численные схемы и программы обеспечивают вычисления собственных функций и их производных по параметру с точностью порядка 0(кр+1), а также собственных значений, их производных и матричных элементов 0(Ь?Р) по шагу Н конечно-элементной сетки [10,11]. Этот факт позволяет использовать вычисленные параметрические собственные значения, собственные функции (базисные функции) и матричные элементы для численного решения с заданной точностью задачи на связанные состояния и задачи рассеяния для различных двумерных эллиптических уравнений, которые сводятся к системам радиальных дифференциальных уравнений в рамках МК [7,19]. Обобщение разработанного алгоритма для решения параметрических двумерных краевых задач в рамках проекционных методов и МКЭ на основе [17,19], которое необходимо при решении краевых задач для трёхмерных эллиптических уравнений, будет дано в следующих работах.
Литература
1. Chuluunbaatar O. et al. Calculation of a Hydrogen Atom Photoionization in a Strong Magnetic Field by using the Angular Oblate Spheroidal Functions // J. Phys. A. — 2007. — Vol. 40. — Pp. 11485-11524.
2. Chuluunbaatar O. et al. Explicit Magnus Expansions for Time-Dependent Schrodinger Equation // J. Phys. A. — 2008. — Vol. 41. — Pp. 295203-1-25.
3. Сечение реакции двух заряженных частиц в канале кристалла / П. М. Кра-совицкий, С. И. Виницкий, А. А. Гусев, О. Чулуунбаатар // Известия РАН. Серия «Физическая». — 2009. — Т. 73. — С. 233-235.
4. Demkov Y. N., Meyer J. D. A Sub-Atomic Microscope, Superfocusing in Channeling and Close Encounter Atomic and Nuclear Reactions // Eur. Phys. J. B. — 2004. — Vol. 42. — Pp. 361-365.
5. Kazaryan E. M., Kostanyan A. A., Sarkisyan H. A. Impurity Optical Absorption in Parabolic Quantum Well // Physica E. — 2005. — Vol. 28. — Pp. 423-430.
6. Чулуунбаатар О. Вариационно-итерационные алгоритмы численного решения задачи на связанные состояния и задачи рассеяния для систем связанных радиальных уравнений // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2008. — № 2. — С. 40-56.
7. Чулуунбаатар О. Многослойные схемы для численного решения нестационарного уравнения Шредингера методом конечных элементов // Вестник РУДН. Серия «Математика. Информатика. Физика». — 2008. — № 3. — С. 68-83.
8. Тихонов А. Н. и Самарский А. А. Уравнения математической физики. — М.: Издательство Московского Унверситета, 1999.
9. Канторович Л. В. и Крылов В. И. Приближенные методы высшего анализа. — М.: Гостехиздат, 1952.
10. Strang G., Fix G. J. An Analysis of the Finite Element Method. — New York: Prentice-Hall. Englewood Cliffs, 1973.
11. Bathe K. J. Finite Element Procedures in Engineering Analysis. — New York: Prentice Hall. Englewood Cliffs, 1982.
12. Abrashkevich A. G., Kaschiejj M. S., Vinitsky S. I. A New Method for Solving an Eigenvalue Problem for a System of Three Coulomb Particles within the Hy-perspherical Adiabatic Representation // J. Comp. Phys. — 2000. — Vol. 163. — Pp. 328-348.
13. Abrashkevich A. G., Abrashkevich D. G., Shapiro M. HSTERM: A program to calculate potential curves and radial matrix elements for two-electron systems within the hyperspherical adiabatic approach // Comput. Phys. Commun. — 1995. — Vol. 90. — Pp. 311-339.
14. Пузынин И. В. и др.. О методах вычислительной физики для исследования моделей сложных физических процессов // ЭЧАЯ. — 2007. — Т. 38. — С. 144232.
15. Gantmacher F. R. The Theory of Matrices. — USA: AMS, Providence, 2000.
16. Голуб Д., Ван Лоун Ч. Матричные вычисления. — М.: Мир, 1999.
17. Chuluunbaatar O. et al. PAROBP: A Program for Computing Eigenvalues and Eigenfunctions and Their First Derivatives with Respect to the Parameter of the Parametric Self-Adjoined Sturm-Liouville Problem // Comput. Phys. Commun. — 2009.
18. Chuluunbaatar O. et al. Benchmark Kantorovich calculations for Three Particles on a Line // J. Phys. B. — 2006. — Vol. 39. — Pp. 243-269.
19. Chuluunbaatar O. et al. KANTBP: A Program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channel Hyperspherical Adiabatic Approach // Comput. Phys. Commun. — 2007. — Vol. 177. — Pp. 649675.
UDC 519.633.2, 517.958, 530.145.6
Algorithm of Numerical Solving the Parametric Sturm-Liouville Problem and Calculation of Solution Derivatives with Respect to the Parameter via the Finite-Element Method
O. Chuluunbaatar
Joint Institute for Nuclear Research Joliot-Curie str. 6, Dubna, Moscow Region, Russia, 141980
The economical algorithm is presented for numerical solving with the given accuracy of the parametrical Sturm-Liouville problem with the third type boundary conditions on the finite interval by the finite-element method with automatical shift of the spectrum. The algorithm includes the calculation with given accuracy a set ~ 10-50 of the eigenvalues, eigenfunctions and their first derivatives with respect to the parameter, and integrals — matrix elements between eigenfunctions and their derivatives with respect to the parameter. Efficiency of the algorithm is demonstrated by numerical analysis of the parametric integrable problem with the parametric third type boundary conditions.
