CC BY
32
Вестник РУДН, Серия Математика. Информатика. Физика. № 2. 2008. с. 49-64 49
УДК 519.633.2,517.958: 530.145.6
Вариационно-итерационные алгоритмы численного решения задачи на связанные состояния и задачи рассеяния для систем связанных радиальных уравнений
О. Чулуунбаатар
Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, Дубна, Московская обл., Россия, 141980
Представлены вариационно-итерационные алгоритмы численного решения задачи на связанные состояния и задачи рассеяния для систем связанных радиальных уравнений в рамках метода Канторовича (МК). Редукция краевых задач с условиями третьего рода для систем связанных радиальных уравнений выполнена методом конечных элементов (МКЭ) высокого порядка точности на неравномерной сетке. В качестве теста для проверки скорости сходимости разложения МК и эффективности аппроксимации задачи МКЭ используются точные значения энергии, фазы и длины рассеяния для модели трёх тождественных частиц (бозонов) на прямой, взаимодействующих парными потенциалами нулевого радиуса. Выполнено сравнение скорости сходимости МК и метода Галёркина в численных расчётах энергии основного состояния данной модели.
Ключевые слова: вариационно-итерационные алгоритмы, метод Канторовича, системы обыкновенных дифференциальных уравнений второго порядка с переменными коэффициентами, метод конечных элементов.
1. Введение
В атомной физике используется гиперсферический адиабатический подход к решению краевых задач для конечномерного уравнения Шрёдингера в области И^ \ {0}, описывающего, например, атом гелия [1] или водородоподобный атом в магнитном поле [2]. В этом подходе, известном в математической литературе как метод Канторовича приведения к обыкновенным дифференциальным уравнениям, уравнение Шрёдингера сводится усреднением по угловым переменным О £ Б^-1 к системе обыкновенных дифференциальных уравнений второго порядка по радиальной переменной р £ И+ с переменными коэффициентами. Переменные коэффициенты выражаются через интегралы между собственными функциями по угловым переменным О вспомогательного оператора угловой части гамильтониана, зависящего от параметра — радиальной переменной р, и их производными по параметру.
Для того, чтобы гарантировать высокую точность численного решения задачи МК, необходимо, прежде всего, иметь алгоритм вычисления матричных элементов с заданной (компьютерной) точностью. Соответствующая программа для вычисления матричных элементов атома водорода в магнитном поле представлена в [2]. Следует отметить, что переменные коэффициенты системы обыкновенных дифференциальных уравнений и её линейно независимые матричные решения при больших значениях радиальной переменой представляются в виде асимптотических медленно сходящихся разложений по обратным степеням по радиальной переменой. Поэтому требуется разработка специальной вариационной постановки задачи, включающей краевые условия третьего рода, и соответствующих стабильных и высокоточных вариационно-итерационных алгоритмов численного решения алгебраических задач, полученных на основе МКЭ.
Статья поступила в редакцию 18 декабря 2007 г.
Работа выполнена в рамках темы ОИЯИ «Математическая поддержка экспериментальных и теоретических исследований, проводимых ОИЯИ 09-6-1060-2005/2009».
Автор благодарит С. И. Виницкого, А. А. Гусева, М. С. Касчиева, И. В. Пузынина и Л. А. Севастьянова за сотрудничество и поддержку.
В данной работе представлены алгоритмы численного решения краевых задач с условиями третьего рода для систем обыкновенных дифференциальных уравнений второго порядка с использованием МКЭ и аппроксимации матричных решений лагранжевыми элементами высокого порядка точности на неравномерной сетке узлов. Обобщённая алгебраическая задача ЛЕ = ЕБЕ относительно пары неизвестных (Е, Е), полученная редукцией задачи на собственные значения для систем дифференциальных уравнений с помощью аппроксимация МКЭ, решается методом обратных итераций в подпространстве с помощью подпрограммы ЯЯРЛСЕ [3]. Обобщённая алгебраической задача (Л — Е Б) Е = Л Ю Е относительно пары неизвестных (Л, Е), полученная редукцией задачи рассеяния в открытых каналах при фиксированных значениях энергии Е, решается с помощью факторизации симметричной матрицы и метода обратной подстановки, используя подпрограммы БЕСОМР и ИЕББЛК [3]. В качестве теста эффективности работы представленного алгоритма и реализующего его программы, а именно: для проверки скорости сходимости разложения МК и эффективности аппроксимации задачи МКЭ используются точные значения энергии, фазы и длины рассеяния для модели трёх тождественных частиц (бозонов) на прямой, взаимодействующих парными потенциалами нулевого радиуса.
2. МК для уравнения Шрёдингера
В ряде случаев конечномерные квантово-механические задачи приводятся к решению стационарного уравнения Шрёдингера для волновой функции Ф(р, О):
(И + и(р,О))Ф (р,О) = ЕФ (р,О), (1)
где И — гамильтониан в ^-мерном пространстве (^ > 1) с эффективным потенциалом и(р, О), Е — спектральный параметр (энергия системы), р = ^Р2 — гиперрадиус системы р £ И/^, О = {Оj}^=1 — набор угловых координат на гиперсфере 5 ^-1(О).
В МК парциальная волновая функция Ф(р, О) ищется в виде разложения по однопараметрическому набору гиперповерхностных функций {Bj (О; р)^=1х:
jmax
Ф(р,О) = £ В,(О;р)х^(р). (2)
j=1
В разложении (2) вектор-функции хМ(р) = (х1г)(р),..., Х(тах(р))Т являются искомыми величинами, а поверхностные функции Б(О; р) = (В1(О; р),...,В.,тах (О; р))Т образуют ортонормированный базис Б(О; р) £ по набору угловых переменных О для каждого значения гиперрадиуса р, который рассматривается здесь как параметр. Базисные функции Bj (О; р) определяются как решения параметрической задачи на собственные значения
(-р2А-Ъ + 2 и(р, О)) В,(О; р) = е,(р)В,(О; р), (3)
где ЛЪ — самосопряжённый оператор обобщённого углового момента. Собственные функции задачи (3) удовлетворяют краевым условиям по угловым переменным О, аналогичным тем, которые наложены на волновую функцию Ф(р, О) при каждом фиксированном значении параметра р £ И,+ :
^Вг(О; р)В, (О; р)^ =| Вг(О; р)В, (О; р)ёО = %, (4)
где 5^- — символ Кронекера, а черта обозначает комплексное сопряжение.
После варьирования подходящего функционала с использованием разложения (2) уравнение (1) сводится к системе из ,]тах обыкновенных дифференциальных уравнений второго порядка относительно неизвестных функций х(р) =
Х(г)(р):
1 „ ё л 1 ё
(L - I) Х(р) S|- 1 + V(p)
+ + Т-Т^ - 2£Ох(Р) = 0. (5)
Здесь I, У(р) и ^(р) — матрицы размерности ^тах х ,?тах, элементы которых определяются соотношениями
^(Р) + (Р)
щ (Р) = (Р) +--^-' нз = '
Vij (р) = Hij (р) +-2-5и , tu = к
= {вЩМ_ | _ Qj (р,_ - (Bi(ß; р, | jM) (6)
При численном решении задачи (1) применяются обобщённые однородные краевые условия третьего рода для парциальной функции Ф'Др, Ю) по гиперрадиусу в граничных точках конечного интервала 0 ^ р < ртах < то:
ß - Л1 ®i(p,ß) = 0, р = 0, ß е Sd-1 (ß), (7)
Р2 ^^^ - Л2 ^(P,ß) = 0, р = Pmax, ß е Sd-1(ß), (8)
др
где ^i, Л1, — вещественные константы, р2 = (ртах) и Л2 = Л2 (ртах) — вещественные числа, зависящие от ртах и определяемые из асимптотических разложений решения. Поскольку гиперповерхностные функции образуют базис, то выполняется условие
(вз(ß; р) ш ß) - Лг ^¿(р, ß)^ = 0, l = 1, 2, (9)
из которого следуют матричные однородные краевые условия третьего рода
№ (i- Qfc)) х(р) - Лгх(р) = 0, l = 1, 2. (10)
Отсюда при l = 1, граничное условие, наложенное на функции х(р) в точке р = 0, может иметь одну из следующих форм:
1) граничное условие Дирихле, если min lim 11Vjj(р)| = то,
ax
Х(0) = 0; (11)
2) граничное условие типа Неймана, если min lim 11 Vjj (р)| < то,
jToр'"1 (ldd - х(р) = 0. (12)
2.1. Граничные условия для задачи дискретного спектра
Для связанных состояний вычисляются энергии Е и радиальная волновая функция х(р)- При больших р радиальная волновая функция х(р) убывает по экспоненциальному или степенному (в случае так называемых полусвязанных состояний) закону. Из уравнения (10) при I = 2 граничное условие для функции х(р) в точке р = ртах может иметь одну из следующих форм:
1) граничное условие Дирихле при = 0
х(р) = 0; (13)
2) граничное условие типа Неймана при Л2 = 0
(IА - Р(р)) х(р) = 0; (14)
3) однородное граничное условие третьего рода при ^2 = 0 и Л2 = 0
(1 А - д(р)) х(р) = л(р)х(р), (15)
где Л(р) = Л2(р)/р2(р). После подстановки (15) в (5) получаем следующую задачу на собственные значения при р = ртах
(у(р) + д2(р))х(р) = Мр)х(р), (16)
где наименьшее собственное значение р(р) и Л(р) удовлетворяет соотношению
1 ёрЛ-1Л(р) рЛ-1 ёр
Заметим, что собственное значение р(р) при р = ртах не зависит от Е и Л(р), что важно для дальнейших вычислений.
= + Л2(р) + 2Е. (17)
2.2. Граничные условия для задачи непрерывного спектра
Для большинства физических задач матричные потенциалы "У(р) и Р(р) при р — ж имеют следующее асимптотическое поведение
V* (р) = Е "рЪ, ^(р) = Е *, для г = 3, (18)
г=2 р г=1 р
V* (р) = * - +^^ + е *, (19)
р р г=з р
где 61 ^ ... ^ — пороговые энергии. В этом случае радиальные волновые функции {х(г) (р)имеют следующее асимптотическое поведение
,) (р) ^ 8шК-(р)) + со8К-(р)) ^ + о /(*+1)М , =1.....N0, (20)
Х*г)(р) - + 0(р-(Л+1)/2 ехр(—(р))) , з = N + 1,...,3тах, (21)
... „Л-1 V /
рЛ-1
где
^(р) = к*р + к* 1п(2к'р) - -п + ¿С,
¿С = arg Г ^--1 к^у) , (р) = Р - ^ 1п(%Р).
(22)
Здесь N — число открытых каналов, ¿с — кулоновские фазовые сдвиги, К = — искомая матрица реакции, к* = у^Е - е* для з = 1,..., N и = л/е* - 2Е для з = N0 + 1, . . . , Зтах.
Теперь рассматриваем квадратичный функционал
ртах
5(Ф, Е, ртах) = I ФТ(р)(Ь - 2Е I) Ф(р)рЛ-1ёр =
0
= П(Ф, Е, ртах) - рт-аХФТ(ртах)Ф(ртах)Л, (23)
где П(Ф, Е, ртах) — симметричный функционал
П (Ф,Е,ртах)= / (^^ ФТ (р)У(р)Ф(р) + 0 ^
+ Фт(р)д(р)ёф|рр) - ^Фр^д(р)Ф(р) - 2ЕФТ(р)Ф(р^ рЛ-1ёр, (24)
и Ф(р) = {ХМ (р)}^ — матрица решения размерностью з'тах х N0, которая удовлетворяет задаче на собственные значения при р = ртах
^ - д(р)Ф(р) = Ф(р)Л, Л = Л«}^. (25)
При использовании МКЭ, уравнение (24) аппроксимируется к задаче на собственные значения при р = ртах (детали см. в разделе 2.4)
с(р)ф(р) = ^ - д(р)Ф(р), (26)
где С(р) — симметричная матрица размерностью з'тах х з'тах. Отсюда получаем соотношение между Ф(р) и её производной ёФ(р)/ёр при р = ртах
^1^ = Я(р)Ф(р), Щр) = С(р) + д(р). (27)
После этого уравнение (27) перепишем с помощью независимых регулярных
Ф^ (р) = {хШр)}^°1 и нерегулярных Ф1ГГ(р) = {^(р)}^ асимптотических решений уравнения (5), а также их производных при р = ртах:
d*(p) = d^reg(p) + d^irr(p)1 dp dp dp
ф(р) = Фreg(p) + Ф*г^) K, ^ = ^^ + K. (28)
Используя формулу (27), приходим к следующему матричному уравнению для матрицы реакций К
(ёФ1гг(р) (ёФгеё(р)
dp
- R(p^i„(p)) K = - (- R(p^reg(p)) . (29)
Регулярные и нерегулярные функции должны удовлетворять обобщённому соотношению при больших p
Wr(Q(p); Ф^), Фreg (p)) = Ioo, (30)
где Wr(•; а(р), Ь(р)) - обобщённый Вронскиан
а(р), Ь(р)) = р
ё-1
аТ(Р) (Т -Ь(Р>) - (^ -а<Р))ТЬ(Р)
(31)
и 100 — единичная матрица размерностью N х При р = ртах соотношение (30) эквивалентно
Wr(R(р); Ф1ГГ(р), Ф^(р)) = Wr(Q(р); «^„(р), (р)). (32)
Таким образом, для матрицы реакции К получаем следующую формулу
К = -X-1 (ртахЩртах), (33)
где
Х(р) = (- R(р)Фiгг(р))o , У(р) = (""¡р^ - ВДФгее(р)^ (34)
— квадратные матрицы размерностью N0 х N0, элементы которых определяются на открытых каналах.
2.3. Регулярные и нерегулярные матричные решения
Мы будем искать регулярные Фгее(р) и нерегулярные Ф;гг(р) матричные решения уравнения (5) с компонентами хге}е (р) = ^ГДр^..., х5Ш!х;(р))Т и х(г)(р) =
(хи (р),..., Хтах¿(р))^, используя следующие асимптотические формы при больших р
Х* (р) = "Ткр^ ¿о+ "Ткр^ ¿о , (35)
Х1гг(р) = СО^КМ^ СЛ'2) . 8!п(щг(р)^ (36)
Х* (р) Укр^ ¿0 р1 + Ук/-1 ¿0 р1 . (36)
Подставляя выражения (18), (19) и (35), (36) в уравнение (5) и приравнивая коэффициенты при 8ш(и>Др)), соз(да;(р)), а также при одинаковых степенях р, получаем рекуррентные соотношения для определения искомых коэффициентов ^1),
(1,2) (¿,1) (¿,2) (0,1) , (0,1) п (0,2) , и с^ , с^ с начальными данными = о^, с^ = 0, с^ = о^,
8^0,2) = 0 следующими из асимптотик (20) [2].
2.4. Формулировка алгебраических задач МКЭ
Для численного решения задачи Штурма-Лиувилля для уравнения (5) с соответствующими граничными условиями из (11), (12) и (13), (14), (15), (25) приближения высокого порядка в МКЭ [3,4] были разработаны в работе [5] и показана их точность, устойчивость и эффективность для широкого набора квантово-механических систем. Высокоточные вычислительные схемы построены из вариационного функционала Рэлея-Ритца для связанных состояний
^шах
ИЗрё-1ф - А рт-Х Е Хортах) > х
шах шах
Яь(х,£,А)= ^ Е [хЯхЬУ"1^ - Арт-Х Е х2(ртах)
Ршах^'
/.ушах ■
Е х|(р)рё_1^р, (37)
¡,^=1 5 = 1
1
и из вариационного функционала Хюльтена для задачи рассеяния
{ршах / \ ( /
/. ,/тах I I ,/тах
I £ [Х(Н - 2Е№|рт-аХ Е х2(ртах)
-1
(38)
А = 1
Здесь
[х(Я - 2Е)хЬ" = - 2£хг(р)хА(р)<^,
(39)
и «' » обозначает производную по р.
Общая идея МКЭ для решения системы обыкновенных дифференциальных уравнений состоит в том, что интервал Л = [0, ртах] делится на подынтервалы Л/ = [Ра_1,Ра] (Л = иП=1 Л/), называемые элементами.
В каждом из подынтервалов Л/ определяются узлы1 ррг и строятся интерполяционные полиномы Лагранжа
фр,г (Р)
Д (Р - Рр,г)
Р
п
РР
(р/,г ^^
С их помощью определяем набор локальных функций N1 (р):
I = 0,
(40)
^ (Р) =
ФР,0 (р), р е Л1,
0, р е Л1,
фР,г (P), ре Л/,
0, р е Л/,
ФР,р(р), ре Л/,
^+1,0^), ре л/+Ъ
0, р е Л/ и л.
фn,p(р), Р е Лп,
0, р е л„,
(41)
I = г + р(з — 1), г = 1,... ,р — 1,
1 = з'р з = 1,...,п -1
1,
I = пр.
Кусочно-дифференцируемые функции {^Р(р)}^=0, Ь = пр, формируют базис в пространстве полиномов порядка р. Далее каждую функцию хДр) аппроксимируем конечной суммой локальных функций ^Р(р):
Х» = £Х^ ^Р(р), х^ = Х^(рР,г).
(42)
=0
Для связанных состояний, после подстановки выражений (42) в вариационный функционал (37) и его минимизации [3,4], получаем обобщённую алгебраическую задачу на собственные значения для определения решений х^ и на наборе локальных функций N (р):
(АР - = 2Е* ВР х
Л
(43)
1 Используется равномерное разбиение каждого из подынтервалов
Н;
Р
Р
Р;-1 +--г, Н; = р; - р;-1, г = 0 ,...,р.
Р
Здесь М — диагональная матрица с нулевыми элементами, кроме последних ^тах элементов, которые равны р^-;. Если используются граничные условия Дирихле или типа Неймана, то А* = 0.
В случае граничного условия третьего рода А* определяется из условия (17)
1 Ёрё-1А* рё-1 Ёр
V = -, '' + (А*)2 + 2Е*, (44)
где V — наименьшее собственное значение задачи (16). Тогда используется следующая итерационная схема для нахождения А = А*, Е = Е* и х = х*
(Ар - А0-1)М)х0-1) = 2Е(0-1) Вр х(5-1),
(А»)2 = V - 0 - А«-1> - 2Е0-1), (45)
ёр ртах
где
А(0)
и х(0) — начальные приближения. Для решения задачи рассеяния при фиксированной энергии Е, после подстановки выражений (42) в вариационный функционал (38) и его минимизации [3,4], получаем обобщённую задачу на собственные значения для определения решений
ФН = ((Х(1))^,..., (х(Мо))^):
= (Ар - 2ЕВр)ф* = МЛ* = (46)
Уравнение (46) удобно переписать в виде
(^ %)(ф*) - рт-аХ (0 0)(%) л = (ф*). (47)
где искомые решения ф* и ф* = Ф(ртах) — матрицы размерностью (Ь^тах -^тах) х N и ^тах х N0, соответственно. Из (47) получаем явное выражение для
ф*
^ а!
ф* = -(саа)-1 саьф* (48)
■ а — аа^ ^аЬфЬ
задачу на собственные значения относительно ф* и Л
сРь - сра(саа)-1сал ф*=рт-хф*л*, (49)
которая с помощью (25) сводится к уравнению
СРь - С^С!)-^) ф* = р^ (Ёфр* - Q(Рmax)Фh) . (50)
Используя (27), получаем соотношение между ф* и её производной:
Ёфр* = R(Рmax)Фh, R(Рmax) = Р^-х (С?ь - С^Щ-1 С^) + Q(Рmax). (51)
Для того чтобы избежать обращения матрицы Саа при вычислении R(рmax), рассмотрим вспомогательное алгебраическое уравнение
( Саа СРь )( Га ) = рё-1 ( 0 )
С ра Сь Г = Ртах I . (52)
и
Это уравнение имеет решение
РР = (ГР )_1рР рР рР = пЛ_Ч РР РР (ГР )_1рР
= (Гаа) ГаЬ* Ь, РЬ = ртах I ГЬЬ - ГЬа(Гаа) Г«Ь
из которого следует, что искомая матрица И,(ртах) равна
Щртах) = (РР)-1 + д(ртах), а решения вычисляются по формулам (48) и (53)
1
р (рр)-1 ф^,
фreg(рmax) + Фirr(рmax)K, (¿)
(53)
(54)
(55)
где К — матрица реакций, определяемая из (33) и Фгее(р) = {х(е^(р)}^=°1; Фirr(р) =
{х(Гг (р)}^==1 — асимптотические решения, определенные в (28).
Таким образом, мы получили выражение для искомых матрицы И,(ртах), которое не содержит собственные значения Л^ и соответствующие собственные вектора Фп. Это позволяет избежать решения алгебраической задачи на собственные значения с матрицей большой размерности и обращения матрицы Ора, что существенно экономит компьютерные ресурсы.
Матрицы жёсткости АР и массы ВР - симметричные и ленточные матрицы, и ВР — положительно определена. Они имеют вид
АР = £ аР, Вр = £ Ь
/=1
/=1
Здесь локальные матрицы ар и Ьр вычисляются по формулам
(а/ Л
+1
/ )'(р)(ФР,г )'(р) + V (р)фР,9 (р)фР,г (р)
-1 5
+ (р) ФР,,(р)(фр,г)'(р) - (ФР,,)'(р)ФР,г(р) Г" р^¿п, (57)
2
Т.,-
2
+1
(ьрг
Т.,-
= / / (р)Фр,г(р)р'_1 "/¿П
^ /
1
2
(56)
где р = р/_1 + 0, 5Т/ (1 + п), 9,г = 0, ...,р, р, V = 1,...,з'тах. Интегралы (57) вычисляются с помощью гауссовской квадратурной формулы
(а/)«
р г 4
= Е V Т4? (ФР,9 )'(рр )(ФР,г )'(рв ) + V (рр )ФР,9 (рр )ФР,г (рр )
р=0 1
г т 2 1 Т •
+ (рр) /(рр)(ФР,г)'(рр) - (ФР,р)'(рр)ФР,г(рр)] ^ }р^_1 у%,(58)
(ьр)«г = Е V фр,, (РР )ФР,Г (РР )р^_1 у %,
р=0
где рр = р/_1 + 0,5Т/(1 + пр), Пр и и>р, д = 0,... ,р — гауссовские узлы и веса.
Для решения поставленной задачи с помощью ЭВМ используется следующая стратегия: зная матричные элементы "У(р) и Р(р), выбираем конечноэлементную сетку (КЭС), потом вычисляем эти матричные элементы в гауссовских узлах и затем вычисляем соответствующие интегралы (57).
Чтобы решить обобщённую задачу на собственные значения (43), мы выбрали метод итераций в подпространстве [3, 4] для решения алгебраических задач на
собственные значения с симметричными ленточными матрицами большой размерности. Для проверки сходимости используемого метода используется соотношение Рэлея для собственных значений и собственных функций. На практике, для достижения машинной точности (2 ■ 10-16), обычно достаточно 10-16 итераций.
Если матрица Ар — А^М в левой части уравнения (43) не является положительно определённой, то данный метод неприменим, однако, используя сдвиг спектра энергии а, матрицу Ар — А^М можно привести к положительно определённой
Ар = Е^ Вр , Ар = Ар — А^М — аВр. (59)
Заметим, что собственные векторы задач (43) и (59) совпадают, и собственные значения отличаются на постоянный сдвиг Ен = Ен + а.
Для решения (16) для граничных условий третьего рода минимальный сдвиг определить сложно, поскольку величины р(ршах) существенно зависят от р = Ртах. Однако можно использовать следующие нижние и верхние оценки величины
р(ршах):
^шах
ИРшах) — < Е ^ ^ = (У(Ршах) + Р^шах^-, « = 1, . . . , Jmax, (60) 5 = М=5
из которых определяется минимальный сдвиг
/ ^шах ч
а = шл? / — Е 0— 1. (61)
\ ...... /
Тогда уравнение (16) принимает вид
(У(ршах) + О? (Ртах) — а1) х(Ршах) = ^(Ртах)х(Ршах), МРтах) = Д(Ртах) + а.
(62)
В этом случае матрица в левой части (62) положительно определена, что позволяет применить метод итераций в подпространстве. Справедливы следующие оценки погрешности [4]
|Е£ — Ет| < С1|Ет|^2р, ||хт(р) — Хт(р) ||0 < С2 |Ет| НР+1, (63)
где Ет, хт(р) — точные решения, Е^, Хт(р) — соответствующие численные решения, Н — максимальный шаг сетки, т — номер решение, С1 и С2 — положительные константы не зависящие от шага Н. Подобные оценки верны также для численных значений (А(г))^ и (х(г)(р))^.
3. Точно решаемая двумерная модель
Краевая задача для модели трёх тождественных частиц (бозонов) на прямой, взаимодействующих парными потенциалами нулевого радиуса, имеет вид [6]:
/1 д д 1 д2 \ Ч ^'д? + Р2 Ир) Ф <Р-9> = 2Е^>- (64)
с граничными условиями на каждом из шести интервалов (вп ^ в < вп+1)
= (—1)к-«с.Ф (р,вк), к = п,п +1,
р дв
Ф(р, в„+1 — 0) = Ф(р, в„+1 + 0), (65)
г дФ (р,в)
11Ш р---= 0.
др
2ЕГас1 = —4с2-2 = —с2* 2Е^х&с1 = —с2К2 = —с2^ . (66)
Здесь вп = к(2п — 1), п = 0 — 5, к = п/6 - параметр, задающий эффективную константу связи д = \f2cR, парного потенциала нулевого радиуса [6]. В случае притяжения (с < 0) у каждой пары системы трёх бозонов имеется парное связанное состояние фо(п) = ^Кехр(—к|п|) с удвоенной энергией — еО0 = с2К2, т.е. 2Е = д2 + е0О), где д — относительный импульс третьей частицы относительно связанной пары [6]. В этом случае известны значения энергий для основного и полусвязанного состояний
¡-2 п2 г 2Еехас1 = с2К2 = с2п
9, ^ = —ск= —сзб
В случае отталкивания (с > 0) ни у одной пары частиц связанного состояния не
существует, е0О) = 0. Для непрерывного спектра в случаях притяжения и отталкивания известна матрица рассеяния Б [6,7]. В случае притяжения фазовый сдвиг
для упругого рассеяния в интервале 0 < д = ^2Е — еО0 ^ —еО0 имеет вид
гехас,(д) = у — а.«ап (). (67)
В силу (67) длина рассеяния равна а0 = — 11ш(1ап(£ехас1(д))/д) = —то. Приведён-
д^О
ные формулы для точных значений энергии, фазы и длины рассеяния используются ниже в качестве теста для проверки скорости сходимости разложения Канторовича и эффективности аппроксимации задачи МКЭ.
3.1. Разложение Канторовича
Рассмотрим формальное разложение решения Ф^ (р, в) уравнения (64), используя набор одномерных базисных функций Вз-(в; р) £ ^21(—п/6, 2п — п/6):
з
шах
Ф (р,в)= Е В (в; р)хЗ (р). (68)
3=1
В разложении (68) функции х(г)(р) = (х1^(р),..., %32ах(р))Т являются неизвестными, а поверхностные функции В(в;р) = (В1(в;р),...,В3шах(в;р))т образуют ортонормированный базис по независимой угловой переменной —п/6 < в ^ 2п — п/6 для каждого значения р, которое рассматривается здесь как параметр. Базисные функции Вз-(в; р) определяются как решения параметрической задачи на собственные значения
1 д2Вз (в; р) р2 дв2 = °з (р)В
ез(р)Вз(в; р), в„ < в<в„+1, п = 0 — 5,
1 = (—1)к-ПскВз(вк; р), к = п, п + 1, (69)
р дв
В (вп+1 — 0; р) = Вз (в„+1 +0; р). Собственные функции Вз-(в; р) удовлетворяют условию нормировки
5 Чп+1
^(в; р)Вз(р; в^ = ]Т I Вг(р; в)В(р; в)ёв = ^. (70)
п=0 л
В результате усреднения уравнения (64) по базису Вз-(в; р) получаем систему обыкновенных дифференциальных уравнений (5) с переменными коэффициентами (6) при ^ = 2 относительно неизвестных функций х(р) = Х(г)(р). Как показано
в [6], краевая задача (69), (70) имеет аналитическое решение. При с < 0 имеем
(буО V р
р) =cosh [6y, (, - nn )], ^=- ( ^ )2,
- ж2) + |ж|
Bj (0; р)
\
у2 + ж2
cos
n(y2 + ж2) - |ж|
6j 0 - ^
, (р)=I j1
2
(71)
где ж = спр/36, целое число n выражается через |0 - nn/3| < п/6, n = 0 - 5. В случае с > 0 используется только формула (71), но индексы j начинаются от 1. Из (69), (70) следуют уравнения для спектра, в случае притяжения (с < 0)
yi (р) tanh(nyi (р)) = -ж, 0 < yi (р) < то,
3 (72)
yj WW^j (р)) = ж, j - 2 <yj (р) <j - 1, j = 2, ...,jmax,
а в случае отталкивания (с > 0)
yj ^W^j (р)) = ж, j - 1 <yj (р) <j - 1, j = 1, . . . , jmax. (73)
Потенциальные матричные элементы "У(р) и Р(р) определяются выражениями (6) и вычисляются аналитически с помощью yj (р) и параметра ж. Их явные выражения и асимптотические разложения (18), (19), а также асимптотические разложения (35), (36) регулярных и нерегулярных матричных решений задачи рассеяния приведены в [8].
3.2. Численные результаты
Изучим сходимость численного решения по МК к известному аналитическому решению по числу уравнений системы (5) для состояний дискретного и непрерывного спектров. В случае притяжения при с = -1 волновые функции Xj (р) связанных состояний имеют асимптотики при больших р [8]
Х1(р) ^ р-1/2ехр (-^р^ Xj(р) ^ Cjр-3ехр (-^р). (74)
Здесь q2 = -2Е + 6i ^ 0, и Cj — константа. Тогда можно использовать граничное условие Дирихле (13) для вычисления энергии основного состояния. Для вычисления энергии полусвязанного состояния q стремится к нулю при близких значениях энергии к точному значению, т.е. функция Х1(р) убывает как р-1/2 и для значений ртах, используемых в численных вычислениях, величину Х1(ртах) нельзя полагать равной нулю. Поэтому здесь используется граничное условие третьего рода (15).
При вычислении энергии основного состояния была выбрана КЭС ^Р[ртш, ртах] = {ртт = 0(1000) ртах = 50} с элементами четвёртого порядка (p = 4). Число в скобках обозначает число конечных элементов на подынтервале. Также было выбрано Зтах = 70 базисных функций. Тогда при Зтах = 70 размерность и ширина ленточных матриц Ap и Bp (56) равны 280040 и 9 соответственно. Расчёты выполнялись на компьютере 2 Alpha 21264, 750 MHz, 2 GB RAM и типом данных с двойной и четверной точностью, обеспечивающим 16 и 32 значащие цифры соответственно.
В табл. 1 показана зависимость погрешности ЛЕь = 2Eh - 2Eexact, численного результата энергий основного состояния 2E|h и её точного значения 2E^xact = -п2/9 от числа уравнений. В первой колонке показано число уравнений Зтах, во второй и четвёртой колонках — численные значения полученных решений на четверной (quad) и двойной (double) точности, в третьей и пятой колонках — погрешность ЛЕ = 2Eh - 2E^xact. Фактор (ж) обозначает 10х. Из табл. 1 видно, что ЛЕь > 0, т. е. получим верхние оценки для точной энергии, что согласуется с
теорией. Кроме того, из табл. 1 видно, что при использование четверной точности МК монотонно с линейной скоростью сходится к точному значению при возрастании ^тах до максимального значения, а при использовании двойной точности — лишь только до 35 уравнений.
Таблица 1
Сходимость численных значений удвоенной энергии основного состояния 2к точному значению 2Е^*3^ = —п2 /9 в зависимости от числа уравнений
jmax Q pquad 2EKM aeTM rt ^double 2EKM л redouble dEKM
1 -1,09644261272635 1, 801(- -04) -1,09644261272582 1, 801(- -04)
4 -1 09662265710171 5,413(- -08) -1,09662265710061 5,413(- -08)
6 -1 09662270528240 5, 949(- -09) -1,09662270528209 5, 950(- -09)
10 -1 09662271083539 3, 967(- -10) -1,09662271083463 3, 979(- -10)
20 -1 09662271122115 1, 099(- -11) -1,09662271122037 1, 245(- -11)
30 -1 09662271123076 1, 390(- 12) -1,09662271122924 3, 276(- -12)
35 -1 09662271123151 6, 357(- 13) -1,09662271122960 3,194(- -12)
40 -1 09662271123182 3, 232(- 13) -1,09662271122940 3, 361(- -12)
50 -1 09662271123204 1, 046(- 13) -1,09662271122827 4,427(- -12)
70 -1,09662271123213 1, 957(- -14)
2 Eexact -1,09662271123215 -1,09662271123215
На примере расчёта энергии основного состояния в табл. 2 представлено сравнение скорости сходимости разложений Канторовича (68) и Галёркина
I I jmax 1
*i(p>0) = Хо(р) + -7= Е Xj(p)cos(6j0). (75)
v 2п vn
Таблица 2
Сравнение скорости сходимости разложений Канторовича (К) (68) и Галёркина (Г) (75) на примере расчёта удвоенной энергии 2Е(]тах) основного состояния
jmax 1 2 3 4 5 6
ЛЕ K ЛЕ r 1, 801(-4) 9, 662(-2) 2, 762(-6) 4,116(-2) 2, 697(-7) 2, 573(-2) 5, 413(-8) 1, 866(-2) 1, 594(-8) 1, 462(-2) 5, 949(-9) 1, 201(-2)
В первой строке показано число уравнений ^тах, во второй и третьей строках показаны разности ЛЕ(^тах) = 2Е(^тах) — 2Еехас1 вычисленного Е(^тах) и точного значения энергии Е|хас1. Медленная скорость сходимости значений энергии ЛЕг(^тах) при увеличении числа ^тах базисных функций (75) объясняется тем фактом, что ослабить достаточные условия сходимости классического метода Галёркина здесь нельзя.
Из формулы (74) получаем граничное условие третьего рода для полусвязанного состояния при р = ртах
lim ^ = 1 + Л %i(pmax),
р >pmax dp \2pmax ) (76)
lim dx(p)- ( 3
+ q) Xj(Pmax). V pmax /
Р^Ршах ёр \Р
В этом случае было выбрано начальное приближение А(0) = —(1/(2ртах) + (?) и КЭС ^Р[рт;п, ртах] = {ртт = 0(500)50(1500) ртах = 500} с элементами четвёртого порядка (р = 4). Длина рассеяния ао = —а = — Иш(Кц/() вычисляется, используя асимптотические разложения регулярных и нерегулярных решений
, Х'1 (р) системы уравнений (5) при больших р с точностью до о(р-5)
х1Г1ее)(р) = ^ + 979776(4(^т5апхб-р4^р-;тах + 1), х11г)(р) =
Х^(р) = 7776(2; - 3)(-1)'+1 (-^ - ) , (77)
хГ(р) = 7776(2; - 3)(-1)'+1 (-^) .
В табл. 4 показана зависимость от числа уравнений разности ЛЕ^ь = 2Е^Ь -2ЕеХас^ численного значения энергии полусвязанного состояния 2Е^Ь и известного точного значения 2Е£Хас1 = -п2/36. В первой колонке показано число уравнений ;'тах, во второй колонке — численные значения полученных решений, в третьей колонке — разница ЛЕ^ь = 2Е^6 - 2ЕеХас1, а в четвёртой колонке — соответствующие собственные значения Л. В последней колонке показана длина рассеяния ао. Фактор (х) обозначает 10х.
Из табл. 4 видно, что с помощью данного метода получена верхняя оценка точной энергии (т.е. ЛЕ^ь > 0), что согласуется с теорией. Заметим, что сходимость по числу уравнений более медленная чем при вычислении основного состояния. Также показана длина рассеяния, соответствующая этой энергии. При увеличении числа уравнений длина рассеяния стремится к -то.
Рассмотрим задачу упругого рассеяния с одним открытым каналом (0 < д ^ - ск). Фазовый сдвиг ¿^ вычислен на КЭС ^р[рт;п, ртах] = {ртт = 0 (1500) ртах = 300/д} с элементами четвёртого порядка (р = 4). В табл. 3 показана сходимость погрешности фазового сдвига, = ¿ехас1 - между численным фазовым сдвигом ¿^ и его точным значением ¿ехас1 в зависимости от числа уравнений и импульса д. С увеличением числа уравнений до ;тах = 25 значения фазового сдвига сходятся к известным аналитическим значениям (67) с точностью 10-3-10-6. Из табл. 3 видно, что > 0, т.е. получим верхние оценки для фазового сдвига, что согласуется с теорией. Из табл. 3 также видно, что при одном и том же числе уравнений значения фазового сдвига получаются тем точнее, чем меньше импульс д.
Таблица 3
Сходимость погрешности фазового сдвига, Л8 = 5ехас* — 5^, между численным фазовым сдвигом 5^ и его точным значением 5ехас* в зависимости от числа уравнений и импульса q.
Фактор (х) обозначает 10х
Ошах q
0, 002 0,100 0,200 0,300 0,400 0, 500 п/6
1 6,180(-1) 2, 972(-2) 3, 946(-2) 5,353(-2) 6,857(-2) 8,513(-2) 8, 930(-2)
4 1, 704(-3) 2, 074(-3) 4,128(-3) 6,188(-3) 8,250(-3) 1,031(-2) 1, 080(-2)
6 4, 019(-4) 1, 285(-3) 2, 566(-3) 3, 848(-3) 5,131(-3) 6,414(-3) 6, 717(-3)
10 8, 213(-5) 7, 299(-4) 1,459(-3) 2,188(-3) 2,917(-3) 3,647(-3) 3, 819(-3)
15 2, 848(-5) 4, 729(-4) 9,462(-4) 1,419(-3) 1,892(-3) 2,365(-3) 2,477(-3)
20 1, 493(-5) 3,470(-4) 6, 954(-4) 1, 044 (-3) 1,391(-3) 1, 739(-3) 1, 821(-3)
25 7, 743(-6) 2, 680(-4) 5,391(-4) 8,105(-4) 1,080(-3) 1,349(-3) 1,412(-3)
Рассмотрим многоканальную задачу упругого рассеяния при 2Е > 0 в случаях притяжения с = -1 и отталкивания с =1. Были рассчитаны матрицы реакции К с ;тах = 6 при 2Е = 0, 0858443221919622 (д = 0, 6) в случае притяжения и при 2Е = 0, 01 (к = = 0,1) в случае отталкивания на КЭС ^р[рт;п, ртах] = {ртп = 0 (2000) 100 (4000) ртах = 2000} и [ртт, ртах] = {ртт = 0 (2000) 100 (4000) 2000 (2000) ртах = 11000} соответственно, с элементами четвёртого порядка (р = 4). Результаты показаны в табл. 5 и 6. Они хорошо согласуются с расчётами работы [6] и приведёнными там аналитическими формулами. Симметрия матрицы реакции К сохраняется с точностью 10-5-10-9. В этих примерах все ;тах каналов открыты, и влияние каждого канала существенно. Поэтому при увеличении числа
Вестник РУДН, Серия Математика. Информатика. Физика. № 2. 2008. с. 49-64 63
уравнений ^тах получаем более точные результаты для подматрицы размерности меньше чем ^тах матрицы реакции К.
Таблица 4
Сходимость численных значений удвоенной энергии полусвязанного состояния к точному
значению = —п2/36 и длины рассеяния ао в зависимости от числа уравнений
^шах 2Е1ъ ЛЕн ъ Л ао
1 4 6 10 15 -0, 274154341 -0, 274155412 -0, 274155416 -0, 274155421 -0, 274155434 1, 336(-6) 2, 657(-7) 2, 615(-7) 2, 563(-7) 2,431(-7) -6, 232(-4) -1,131(-3) -1,131(-3) -1,129(-3) -1,122(-3) -699,166 -301369, 617 -1338701, 715 -7577671, 004 -27997011, 620
о глеха^ 2ЕНЪ -0, 274155667 -те
Таблица 5
Элементы матрицы реакции К в случае притяжения с = —1 и д = 0, 6. Фактор (х)
обозначает 10х
Л1 1 2 3 4 5 6
1 -0,1261 0, 5000(-7) 0, 9473(-8) 0, 6237(-7) 0, 2544(-6) 0, 6469(-6)
2 0,4885(-7) 0, 5940 0, 3295(-1) 0, 6433(-2) 0, 2276(-2) 0,1024(-2)
3 0, 7318(-9) 0, 3295(-1) 0, 6333 0, 4163(-1) 0, 9727(-2) 0, 3783(-2)
4 0,1269(-9) 0, 6432(-2) 0,4163(-1) 0, 6349 0,4301(-1) 0,1031(-1)
5 0, 6478(-8) 0, 2275(-2) 0, 9726(-2) 0, 4299(-1) 0, 6250 0, 4228(-1)
6 0,1635(-7) 0, 1023(-2) 0, 3781(-2) 0,1030(-1) 0,4221(-1) 0, 5836
Таблица 6
Элементы матрицы реакции К в случае отталкивания с =1 и к = 0,1. Фактор (х)
обозначает 10х
Л ^ 1 2 3 4 5 6
1 —2, 56205 —0, 6155(-1) —0,1278(-1) — 0,4622(-2) —0, 2187(-2) —0,1218(-2)
2 —0, 6155(-1) —2, 6364 —0, 7897(-1) — 0,1960(-1) —0, 8056(-2) —0,4189(-2)
3 —0,1278(-1) — 0, 7897(-1) —2, 6446 — 0,8259(-1) —0, 2175(-1) —0, 9458(-2)
4 —0,4622(-2) — 0,1960(-1) —0, 8259(-1) — 2,6558 —0, 8501(-1) —0, 2319(-1)
5 —0, 2187(-2) —0, 8056(-2) —0, 2175(-1) — 0,8500(-1) —2, 6937 —0, 8912(-1)
6 —0,1218(-2) —0, 4189(-2) —0, 9459(-2) — 0,2319(-1) —0, 8909(-1) —2, 80482
4. Заключение и обсуждение результатов
Представлена новая вариационная формулировка задачи на связанные состояния и задачи рассеяния для системы связанных радиальных уравнений в МК, включающая краевые условия третьего рода. Разработаны стабильные и высокоточные вариационно-итерационные алгоритмы численного решения краевых задач на основе МКЭ. Выполнены численные исследования скорости сходимости к точному решению для точно решаемой модели трёх тождественных частиц на прямой, взаимодействующих парными потенциалами нулевого радиуса. Показано, что при использовании разложения Галёркина (75) медленная сходимость значений энергии ЛЕг (^тах) при увеличении числа ^тах базисных функций объясняется тем фактом, что ослабить достаточные условия сходимости классического метода Галёркина здесь нельзя. Более того, редуцированная таким образом задача при любом конечном ^тах на полуоси имеет ложный кулоновский спектр с квантовым дефектом и качественно отличное пороговое поведение для фазового сдвига. Это обстоятельство следует учитывать при редукции методом Галёрки-на краевой задачи для системы трёх частиц с парными взаимодействиями при наличии связанных состояний.
ф
ф
ф
ф
Литература
1. Abrashkevich A. G., Kaschiev 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. — Vol. 163. — 2000. — Pp. 328-348.
2. Chuluunbaatar O. et al. POTHMF: A Program for Computing Potential Curves and Matrix Elements of the Coupled Adiabatic Radial Equations for a Hydrogen-Like Atom in a Homogeneous Magnetic Field // Comput. Phys. Commun. — 2007.
3. Bathe K. J. Finite Element Procedures in Engineering Analysis. — New-York: En-glewood Cliffs, Prentice Hall, 1982.
4. Strang G., Fix G. J. An Analysis of the Finite Element Method. — New-York: Englewood Cliffs, Prentice Hall, 1973.
5. Finite-Element Solution of the Coupled-Channels Schrödinger Equation Using High-Order Accuracy Approximations / A. G. Abrashkevich, D. G. Abrashkevich, M. S. Kaschiev, I. V. Puzynin // Computer Physics Communications. — Vol. 85. — 1995. — Pp. 40-65.
6. Chuluunbaatar O. et al. Three Identical Particles on a Line: Comparison of Some Exact and Approximate Calculations // J. Phys. A. — Vol. 35. — 2002. — Pp. L513-
7. Kuperin Y. A. et al. Connections and Effective S-matrix in Triangle Representation for Quantum Scattering // Annals of Physics. — Vol. 205. — 1991. — Pp. 330-361.
8. Chuluunbaatar O. et al. KANTBP: A program for Computing Energy Levels, Reaction Matrix and Radial Wave Functions in the Coupled-Channels Hyper-Spherical Adiabatic Approach // Comput. Phys. Commun. — Vol. 177. — 2007. — Pp. 649675.
UDC 519.633.2,517.958: 530.145.6
Variational-Iteration Algorithms of Numerical Solving Bound State and Scattering Problems for Coupled-Channels Radial
Equations
The variational-iteration algorithms of the numerical solving of the bound state and scattering problems for coupled-channels radial equations are presented in the framework of the Kantorovich method (KM). Reduction of the boundary problems with conditions of the third type for systems of coupled radial equations is executed by a finite element method (FEM) with a high order accuracy on a non-uniform grid. As benchmark calculation to check rate of convergence and decomposition KM and efficiency of the FEM approximations of problem, we used exact values of energy, a phase and lengths of dispersion for model of three identical particles (bosons) on a straight line interacting by the pair zero-range potentials. A comparison of convergence rate of KM and Galerkin method in numerical calculations of the ground state energy of the given model is performed.
L525.
O. Chuluunbaatar
Joint Institute for Nuclear Research Joliot-Curie str., 6., Dubna, Moscow Region, Russia, 141980
Ф-
Ф
-Ф
Ф
