CC BY
50
Вестник РУДН, Серия Математика. Информатика. Физика. № 1. 2008. с. 43-53 43
УДК 519.633.2, 517.958: 530.145.6
Многослойные схемы для численного решения нестационарного уравнения Шрёдингера
О. Чулуунбаатар
Объединённый институт ядерных исследований Россия, 141980, Дубна, Московской обл., ул. Жолио-Кюри, 6
Представлен алгоритм генерации в среде MAPLE и REDUCE многослойных неявных схем для численного решения нестационарного уравнения Шрёдингера на основе факторизации унитарного оператора эволюции. Исследуются оптимальные методы построения дополнительных калибровочных преобразований, позволяющие выделять симметричные операторы, необходимые для генерации экономичных алгебраических эволюционных схем по пространственной переменной методом конечных элементов. Эффективность сгенерированных схем до шестого порядка точности по шагу временной переменной и до седьмого порядка точности по шагу пространственной переменной демонстрируется численным анализом интегрируемой модели осциллятора во внешнем поле, зависящем от времени.
Ключевые слова: нестационарное уравнение Шрёдингера, задача Коши, разложение оператора эволюции, калибровочные преобразования, операторно-разностные многослойные схемы, метод конечных элементов.
Введение
С развитием фемтосекундной и тераваттной лазерной физики компьютерное моделирование динамики атомных систем в поле мощного лазерного импульса требует развития новых схем решения задачи Коши на конечном интервале для нестационарного уравнения Шрёдингера [1]. Кроме того, развитие нанотехноло-гий привлекает внимание к решению задач управления квантовыми системами, которые описываются нестационарным уравнением Шрёдингера или соответствующим квантовым уравнением Лиувилля [2]. Новые схемы необходимо также развивать для решения эволюционных задач, описывающих в рамках метода функционала плотности для многоэлектронных систем в поле ядер и в дополнительных электромагнитных полях [3], как тормозное излучение при а—распаде тяжёлых ядер [4], так и процессы каналирования ионов в кристаллических плёнках [5,6].
В этом случае возникает потребность разработки эффективных алгоритмов и многослойных схем [7,8], учитывающих специфику задачи в рамках методов расщепления [9] и дискретных схем по пространственной переменной [10-12]. Для эффективного применения метода конечных элементов необходимо исследовать дополнительное калибровочное преобразование эффективных гамильтонианов многослойных схем высокого порядка точности по шагу временной переменной, которые для вещественных гамильтонианов редуцируются к алгебраическим задачам с симметричными матрицами на локальных носителях.
В работе построены и исследованы многослойные схемы для численного решения нестационарного уравнения Шрёдингера до шестого порядка точности по временной переменной и до седьмого порядка точности по пространственной переменной. Для факторизации оператора эволюции и генерации неявных схем решения нестационарного уравнения Шрёдингера с заданной точностью как
Работа выполнена в рамках темы ОИЯИ «Математическая поддержка экспериментальных и теоретических исследований, проводимых ОИЯИ 09-6-1060-2005/2009'».
Автор благодарит С. И. Виницкого, В. П. Гердта, А. А. Гусева, В. Л. Дербова, М. С. Касчиева, И. В. Пузынина, В. А. Ростовцева и Л. А. Севастьянова за сотрудничество и поддержку. Статья поступила в редакцию 15 октября 2005 г.
по временной, так и пространственной переменным применён символьный алгоритм [13], реализованный в среде MAPLE и REDUCE. Структурные элементы алгоритма реализованы на основе методики, разработанной в [14,15]. Соответствующие эффективные гамильтонианы и алгебраические матричные задачи генерируются на основе методов расщепления для оператора эволюции на равномерной сетке по времени [16-18] и конечноэлементных сетках по пространственной переменной в рамках МКЭ [11,12], используя подходящие интерполяционные полиномы [19,20].
Содержание работы следующее. В разделе 1 представлен алгоритм вычисления оператора эволюции для решения задачи Коши нестационарного уравнения Шрёдингера. В разделе 2, используя Паде-аппроксимацию оператора эволюции, приводится алгоритм генерации многослойных разностных схем решения задачи по временной переменной. В разделе 3 изучается дополнительное калибровочное преобразование эффективных гамильтонианов многослойной схемы, необходимое для эффективного применения метода конечных элементов по пространственной переменной. В разделе 4 метод конечных элементов используется для конструирования численных схем требуемого порядка точности относительно шага конеч-ноэлементной сетки по пространственной переменной. В разделе 5 стабильность и эффективность сгенерированных алгебраических схем показана на примере расчёта модельной одномерной задачи. В заключении обсуждаются основные результаты и дальнейшие перспективы применения разработанного алгоритма и сгенерированных многослойных схем.
1. Вычисление оператора эволюции
Рассмотрим нестационарное уравнение Шрёдингера с гамильтонианом H(ж, t), описывающее динамику атома во внешнем поле, на интервале времени t g [to,T] для начального состояния ^o(x):
д^(ж, t)
г-—— = H(x,t)^(x,t), ^(x,to) = ^o(x),
dt
= J |^(ж, t)|2dx = 1,
(1)
^(x,t) g Wl(Rn) ® [to,T], ^o(x) g Wl(Rn).
Для решения задачи используется равномерная сетка по временной переменной £
ПТ [£о,Т ] = {£о, = ¿к + = Т} (2)
с шагом т = £к+1 — ¿к, (к = 0,1,.. .К). Решение в момент времени £ свя-
зано с решением ^(ж,^) в момент времени ¿к унитарным оператором эволюции
и(Мк, А):
.ди,А) = АЯ(х,^)^(£,£к, А), и(£к,4,А) = 1,
где А — формальный параметр. Решение задачи ^(¿к+О = ^(ж^к+О на к + 1-м слое определяется через решение ^(¿к) на предыдущем к-м слое с помощью оператора эволюции и(¿к+1, ¿к, А):
^(¿к+1) = и(£к+1, ¿к,А)^(£к), и(£к+1, ¿к, А) = ехр{—ггАк^к+ъ А)}. (3) Здесь операторы Ак (¿, А) представляются в виде ряда
Ак(¿, А) = - £ А*A(J)k(¿), Ак(¿к, А) = 0, (4)
т *=1
по степеням формального параметра А. Неизвестные операторы А^)^) = А(*)к(¿) вычисляются с помощью алгоритма, представленного ниже, после чего формальный параметр А полагается равным единице (А = 1).
2
1. Коэффициенты ряда (4) вычисляются из тождества [21]
-■ан(ь) = ^^ ]Г + (айа^ь))... (айА(1ч)(ь))А(П)(ь). (5)
п= 1 9=0 11,...,1д = 1 + )!
Здесь линейный оператор (айА) : С(Х) ^ С(Х) определён для операторов Л, В £ С(Х) в виде (айА)В = [А, В ] = АВ — В А и имеет следующие свойства: (айА)0 В = В, (айА)В = (айА)*-1 (айА)В; точкой над оператором А(п)(Ь) обозначена частная производная А(п)(Ь) = д4А(п)(Ь) по временной переменной в момент времени Ь.
2. Приравнивая коэффициенты при одинаковых степенях параметра У в левой и правой частях тождества (5), получаем рекуррентную систему дифференциальных уравнений первого порядка
■н (ь) + А {1)(г) = 0, 2(айА(1)(г))А (1)(г) + А (2)(ь) = 0,
2(айА(2)(г))А (1)ф + 6(айА(1)(г))2А (1)(г)+ (6)
+ 1(айА(1)(г))А (2) (г) + А (3)(ь) = 0,
3. При заданном Н (Ь) из системы (6) находим эффективные операторы А0)(Ь) = А(*)к(Ь):
tk + 1
1
А{1)(г) = - I ¿ьн(ь),
tk+1 4
А(2)(г) = ^ J М!(айн(ь''))н(ь')
4к+1 4 4' (7)
А{3)(Ь) = — ¿Ь(айН(Ь))I ¿Ь'(айН(Ь'))! ¿Ь''Н(Ь'')
4к 4к 4к
— ^"¡¿ь (ай ¿ь'н (Ь')^ Н (ь), 4к V V" ) )
Отметим, что операторы А(*)(Ь) связаны с известным гамильтонианом Н(Ь), зависящим от временной переменной Ь, разложением Магнуса [22].
4. В предположении, что оператор Н(Ь) = Н(х,Ь) зависит гладко от Ь и от х, с помощью разложения
1 2М
Акм) = А(кМ)(Ьк+1) = - £ А*А(Лк(Ьк+1), (8)
т *=1
получаем аппроксимацию унитарного оператора эволюции (3) с точностью порядка 0(т2М+1) относительно шага т сетки (2),
и (Ьк+1,Ьк; т) = и(м )(Ьк+1,Ьк; т) + 0(т2М+1) = ехр(—гтАкМ)) + 0(т2М+1). (9)
е—
е
—е е
46 Чулуунбаатар О. Многослойные схемы для численного решения. . .
5. Разлагая гамильтониан Н(¿) в ряд Тейлора в окрестности точки ^ = ¿к + т/2
Н(¿) = Е ( — ^1/2)* д*Н(¿к+1 /2), (10)
вычисляя интегралы в (7) и полагая формальный параметр А равным единице (А = 1), получаем выражения для операторов Ак"\ При М = 1, 2, 3 операторы АкМ\ аппроксимирующие оператор эволюции согласно (8) и (9), имеют вид:
А<м > = /¡к" > + 1.4<м >,
А1 А1
А?
A3
= H, = 0,
= Â(i) + T2 H
= A +24 H
= + 12 (adH )H,
(11)
4 4 4
+ I92ÔH — 72o(adH )2H - ¿щ^»^
4 ... 4 4
aA ï3) = A ï2) - 48ô(adH.)H + isô^)^ + Тй^^,
Здесь гамильтониан H = H(x,ifc+i/2) и его частные производные H = dtH(x,t)|t=tfe 2 по временной переменной t вычисляются в момент времени t = ti+i/2- Заметим, что при разложении гамильтониана (10) в ряд Тейлора в центральной точке интервала t g [ti, ti +i], оператор эволюции не содержит нечётных степеней т, поскольку операторы U(ti,ti+i,A) и U(ti+i,tfc,A), определённые на одном временном интервале, являются взаимообратными.
2. Генерация операторно-разностной схемы
Для генерации из (3) неявной операторно-разностной схемы представим экспоненциальный оператор (9), (11) на каждом k-м слое в виде обобщённого [M/M] Паде-разложения с той же точностью
M т + a(M ),(M )
exp(-iT,<M)) = П + O(T2M+■). Tzï = f + 2M%),?M). (12)
Z=i 1 + 2M Л i
Здесь коэффициенты ), Z =!,•••, M, M ^ 1 являются корнями полиномиального уравнения iFi( —M, —2M,2Mi/a) = 0, где iF вырожденная гипергеометрическая функция, коэффициенты ) = Re aZ-M) — iIm ) комплексно сопряжённые к ). Коэффициенты ) (см. табл. 1) имеют следующие свойства: Ima^M) < 0 и 0, 6 < )| < где ^ ~ 0,28 корень трансцендентного урав-
нения ^exp(^ + 1) = 1. Заметим, что условие т < 2M^||,iM) ||-1 гарантирует
справедливость приближения для всякого ограниченного оператора A( M) .
Используя разложение (12), произведём факторизацию эволюционного оператора, т. е. представим (3) в виде системы M уравнений относительно набора M — 1
вспомогательных функций ф^ :
^fc+C/M = Tf i ф fc+(C-1)/M, Z =1,...,M, (13)
Вещественные и мнимые части коэффициентов ^, М = 1, 2, 3, £ = 1,. .. , М
Таблица 1
м С И,е а 1т а
1 1 0 -1.0
2 1 -0.57735026918962576450914878050 -1.0
2 2 +0.57735026918962576450914878050 -1.0
3 1 -0.81479955424892281841473623156 -0.85405673065166346526579940886
3 2 +0.0 -1.29188653869667306946840118228
3 3 +0.81479955424892281841473623156 -0.85405673065166346526579940886
решение которой обеспечивает переход от ф(1к) на к-м слое к ф(1к+1) на к + 1-м слое. Поскольку Ъпа^^ < 0 вспомогательные операторы Т^к изометрические, следовательно норма всех \\фк+^/Му равна единице с точностью 0(т2М)
\\фк\\ = \\фк+1/М у = ... = \\фк+1\\. (14)
Применение соотношений (13) для приближённого решения задачи Коши (1) на каждом к-м слое сетки (2) позволяет генерировать неявную операторно-разностную схему [16]:
фк = ф(х,1к),
I}4М})й = (I}АкМ))фГ, С = 1,2,...,М, (15)
2Ма Ак )фк- 2Ма ,
Ф(х,гк+1) = ФМ, к = 0,1,...,к —1.
Вспомогательные функции фк, (£ = 1,...,М — 1) в (15) можно рассматривать как некоторые приближённые решения на дробных шагах 1к+^/М = Ък + т^/М, £ = 1,...,М — 1 временного интервала \Ьк,1к+1], соответствующего к-му шагу т сетки (2). Заметим, что неявная унитарная М-слойная схема (15) сохраняет норму (14) разностного решения и является стабильной [10]. Каждое уравнение по отдельности из системы М уравнений (15) имеет порядок аппроксимации не выше второго, тогда как вся схема (15) имеет суммарную аппроксимацию 2М-того порядка (0(т2М)) по шагу сетки (2) в смысле [10]. При 2М = 2 схема (15) сводится диагональной Паде-аппроксимацией ранга [1/1] оператора эволюции к известной схеме Кранка-Николсона [23].
3. Схемы с симметричными операторами
Рассмотрим задачу Коши для нестационарного уравнения Шрёдингера, описывающего динамику одномерного модельного атома с потенциальной энергией V(х) во внешнем поле f (1) в дипольном приближении на интервале 1 £ ]
д 1 д2 1 -ф(х, 1) = Н (х, 1)ф(х, 1), Н (х, 1) = — ^ ^ + f (1)х + V (х),
Ф(—(х>,г) = ф((х>,г) = 0, ф(х,г0) = ф0(х), (16)
сю ■ '
2 _ / „/,*
= I ф*(х,1)ф(х,1)^х =1.
В этом случае операторы (11) имеют вид
А'м > = Лкм 1 + 1-4'" >
.¡к1 4?
-¡к3
-4к3
-¡к41 = -¡к3
1 д2
— 1 дХ2 + /ж + У (Х),
= 0, =.¡к1 =.«к1
= .¡к2
= -«к2
т2 + 24
_т!, ±
12 ^ дж,
+ ж?. — ^//• — ^ / + ^ ¿2
1920 1 720^ 720 дж ^ 240-7 ' т4 ...д2У (ж) • д3У (ж) •
— 480 ^ дж — 720 дж2 ^ дж + 1440 дж3 ^ т8 д3У •• т8 д4У •• д т8 д5У
^ а™2 ополп Я™4 ^ Я™ 1 оппсп а~5 ^
30240 дж3 дж2 30240 дж4 ^ дж 120960 дж5
т8 д2УдУ •• т8 д2V •• т8 д2У -2 т8 дУ.... ---/---/ / +---/2---f +
30240 дж2 дж 30240 дж2 ^ 7560 дж2 ^ 40320 дж ^
т8 ............т8 ________т8 ... • т8
+__f ж__ff +__ff__р
322560 40320 5040^ 24192'' '
Х(4) _ *(3) _ т8 д4У • т8 д5У • _ т8 д6У •
к = к — 30240 дЖ4"/&Г3 — 20160 аХ5^ — 40320 дЖ6^дж —
_ т8 д3У дУ • _ т8 д3У • _ т8 /д^УV / — 10080 дЖ3 аж^дж- 10080 дЁ3"— 30240 VдЖ2/ —
т8 д2У ... д т8 д4У • т8 д3Уд2У-
/ а~ оп 1£п а_4 // юпп<? а~3 а_2 /
40320 дж2 дж 20160 дж4 12096 дж3 дж2 _ т8 д4У дУ т8 ..... _ т8 д7У т8 д3У...
— 20160дЖ4"дж— 53760 ^ дж-241920аХ7^— 80640 дЖ3^
(17)
Здесь / = /(¿к+1 /2), / = (¿к+1 /2),.... Для получения схем с частичным расщеплением оператора эволюции применяем калибровочное преобразование
ф = ехр (^)) ф, Д"1 = ехр (^)) 1 ехр (—)) . (18)
Оператор ) ищем в виде ряда по степеням т
2м
) = Е т 3 (19)
з=о
с неизвестными коэффициентами , которые определяются из условия
¿к") Ф = 0(т2М), Iк")= ехр (1^м)) ¿к") ехр (—)) . (20)
Подставляя разложение (19) в уравнение (20) и приравнивая коэффициенты при одинаковых степенях т, получаем алгебраическую систему рекуррентных уравнений для нахождения искомых коэффициентов с начальным условием £о = 0.
Таким образом вместо (15) получаем неявную операторно-разностную схему
фк = ехр )) ф(х,1к), I )АкМА ф< = (I )АкМА фк-1, С = 1,2,..., М, (21)
2Ма Ак )фк- ТМа
ф(х,1к+1) = ехр (-1БкМ)) фМ, к = 0,1,...,К — 1,
которая использовалась при численных расчётах (см. параграф 5). Соответствующие эффективные операторы Ак^) и БкМ) имеют вид
2
= —1 ^ + fx + V (х) 2
Ак2) = Ак1) + ^ х!,
(3) = а(2) хТ 12 + ^(х) г
' к 1920 1 1440-' 1 14401 720 дх 1, (22)
Б к!1 = 0,
2
Б(2 = Бк1) + Т- х/,
Б(3) = Б(2) _ II ^(х) ' + IIх ...
Бк = Бк 720 дх 1 + 480х'.
Заметим, что ) при М ^ 3 представляет собой оператор умножения. Однако при М = 4 не удаётся избавиться от несимметричного оператора из (17), содер-
- д3
жащего третью частную производную по пространственной переменной д^з, если Бк4) определён как оператор умножения. Поэтому при построении схем более высоких порядков М ^ 4, БкМ) представляет собой оператор общего вида. В этом случае для операторной факторизации экспоненциальных операторов в схеме (21) необходимо использовать формулу Хаусдорфа
с
вхр(А)Б вхр(—А) = ^ -^айА)Б.
3=0
3!
Другая возможность состоит в прямой факторизации (18) исходного оператора эволюции
и (Ьк+1,1к ,Х) = ехр(—гтБк (X)) вхр(—гтСк (X)). 4. Генерация алгебраической задачи МКЭ
Для численного решения задачи (1) на сетке От граничные условия и условие нормировки по пространственной переменной х на бесконечном интервале заменяются соответствующими условиями на конечном интервале хт1п < х < хтах (|хШт| ^ 1, |хШах| ^ 1). Далее, на каждом к-м шаге сетки От мы рассмотрим дискретное представление решений ф(х; 1к) задачи (16) с помощью метода конечных элементов на сетке
^р — [х0 — xm\n, х3 — х3-1 + , хп — хтах], 3 — 1, . . . , 71
в виде конечной суммы локальных функций
пр
ф(х; 1к) = ^ ф„(1к)%(х), (23)
^=0
где фД^ ) являются узловыми значениями неизвестной функции ф(ж; tk ), относительно которых численно решается исходная задача. Локальные функции являются кусочно-непрерывными полиномами порядка p, равными единице только в узлах ж^ и равны нулю в остальных узлах сетки i?p, т.е., N^(xv) = , v = 0,..., np. Коэффициенты ^(tfc) формально связаны со значениями вектора задачи ф(ж; tk) в Лагранжевых узловых точках [24]:
h ■
^M(tk ) = ф(ж[м]; tk ), жм = X[r+p(j-i)| = Xj-1 + -1 r, r = 0,1,...,p.
p
Перепишем каждое из дифференциальных уравнений системы (21) на конечном интервале жшш < ж < жшах с граничными условиями на концах интервала Ф(Жш1п, t) = Ф(Жтах, t) = 0 в виде
A^k = B^-1. (24)
Подставляя разложение (23) в схему (24), умножая слева на N^(x) и интегрируя на интервале А = u"=i^j, получаем систему линейных алгебраических уравнений
(4 )v = V (4-1)v. (25)
Как известно [11], теоретическая оценка по норме || ■ || на сетке i?p между точным ф(ж; tk) и приближенным (tk) решениями имеет порядок (tk) — ф(ж, tk)|| = O(hp+1). Для решения задачи (25) с комплексной симметричной 2p+1-диагональной матрицей использовалась процедура F07QSF (ZSPTRS) из библиотеки Фортрановских программ NAG [25].
5. Модельный расчёт
Задача (16) решалась на интервале 0 ^ t ^ 10 для осциллятора V(ж) = и2ж2/2 с частотой и во внешнем периодическом поле f (t) = fo sin(uot) напряжённостью fo и угловой частотой uo. В качестве начального состояния фо(ж) был выбран гауссов волновой пакет
фо(ж) = 4 u/п exp(—и (ж — жо)2/2 + ipo(ж — жо)).
Были выбраны следующие значения параметров жo = 0, po = 1, u = uo = 1, f = 0,25 и конечно-элементная сетка с 1000 элементами шестого порядка p = 6 на интервале [жшш = —20, жшах = 20], обеспечивающими точность (O(h7)) при постоянном шаге h = 0.04. Вычисления интегралов проводились по квадратурным формулам Гаусса для полиномов Лежандра порядка p +1 = 7. Расчёты выполнялись на компьютере 2 Alpha 21264, 750 MHz, 2 GB RAM, используя Compaq Fortran 77 и тип данных real*16 (complex*32), обеспечивающий 33 значащие цифры.
Для анализа сходимости схем (21)-(22), применяемых для решения задачи, на последовательности трёх вдвое сгущающихся временных сеток вычислялись абсолютная погрешность
xmax
Er2(t,j)= У |ф(ж,t) — (ж^)|2ёж (26)
xmin
и коэффициент Рунге [26]
a(t) = log
Er(t, 1) — Er(t, 2)
Er(t, 2) — Er(t, 3)
(27)
Вестник РУДН, Серия Математика. Информатика. Физика. № 1. 2008. с. 43-53 51
где индекс ] = 1, 2, 3 нумерует решения, полученные для шагов т, т/2, т/4 по временной переменной, а ф(х,1) — известное аналитическое решение [2].
В качестве примера на рис. 1 показаны графики = 1, 2,3 для схем с
М = 1, 2, 3. Верхние три кривые соответствуют схеме (21)-(22) с М = 1, средние — с М = 2 и нижние — с М = 3. Возле каждой тройки кривых показано среднее значение а по всем значениям а(1к) на сетке От [0,10]. В результате подтверждены теоретические оценки порядка сходимости предложенных схем (21)-(22), соответственно, второго (М = 1), четвёртого (М = 2) и шестого (М = 3) порядка точности 0(т2М).
Рис. 1. Логарифм абсолютной погрешности Ег(£; ]), ] = 1, 2, 3 (штрих-пунктирная, пунктирная и сплошная кривые), соответственно для схем второго (2М = 2), четвёртого (2М = 4) и шестого (2М = 6) порядков.
Заключение
Методика, реализованная в виде символьных алгоритмов в среде MAPLE и REDUCE, позволяет генерировать многослойные схемы и соответствующие подпрограммы в среде FORTRAN численного решения эволюционных задач для нестационарного уравнения Шрёдингера на конечном интервале с заданной точностью как по временной, так и по пространственной переменным.
Изучена методика построения дополнительных калибровочных преобразований для построения более экономичных операторно-разностных схем, которые для вещественных гамильтонианов редуцируются к алгебраическим задачам с симметричными матрицами на локальных носителях, что обеспечивает эффективность применения МКЭ.
На модельном расчёте явно показано строгое соответствие численных оценок сходимости сгенерированных схем второго, четвёртого и шестого порядка точности относительно шага равномерной сетки по временной переменной теоретическим оценкам.
Развитая техника может быть обобщена при решении эволюционных задач для нелинейных уравнений Шрёдингера [27,28] с использованием симметрийного формализма групп Ли [29] и итерационных ньютоновских схем [30]. Возможна также адаптация метода для расчёта динамики атомных систем во внешнем лазерном импульсном поле, т.е. для гамильтониана системы, содержащего быстро осциллирующие функции от временной переменной.
ф
ф
0
0
Литература
1. Coulomb-Volkov approach of atom ionization by intense and ultrashort laser pulses / G. Duchateau, E. Cormier, H. Bachau, R. Gayet // Phys. Rev. A. — Vol. 63. — 2001. — Pp. 053411-1-53411-11.
2. Butkovskiy A. G., Samoilenko Y. I. Control of Quantum-Mechanical Processes and Systems. — Dordrecht Hardbound: Kluwer Academic Publishers, 1990.
3. Киржниц Д. А. Полевые методы теории многих частиц. — М.: Госатомиздат,
4. Misicu S., Rizea M., Greiner W. Emission of electromagnetic radiation in a-decay // J. Phys. G: Nucl. Part. Phys. — Vol. 27. — 2001. — Pp. 993-1003.
5. Khodyrev V. A. The treatment of energy loss of induced current density // Advances in quantum chemistry. — Vol. 45. — 2004. — Pp. 124-158.
6. Demkov Y. N., Meyer J. D. A sub-atomic microscope, superfocusing in channeling and close encounter atomic and nuclear reactions // Eur. Phys. J. B. — Vol. 42. — 2004. — Pp. 361-365.
7. Bankov N. G., Kaschiev M. S., Vinitsky S. I. Adaptive method for solving the time-dependent Schrodinger equation // Comptes rendus de l'Akademie bulgare des Scienses. — Vol. 55. — 2002. — Pp. 25-30.
8. Adaptive numerical methods for time-dependent Schroedinger equation in atomic and laser physics / V. L. Derbov, M. S. Kaschiev, V. V. Serov et al // Proc. SPIE. — Vol. 5069. — 2003. — Pp. 52-60.
9. Marchuk G. I. Partial Differential Equations. II SYNSPADE-1970. — New-York, London: Academic Press, 1971.
10. Самарский А. А. Теория разностных схем. — М.: Наука, 1977.
11. Стренг Г., Фикс Д. Теория метода конечных элементов. — М.: Мир, 1977.
12. Bathe K. J. Finite Element Procedures in Engineering Analysis. — New-York: Englewood Cliffs, Prentice Hall, 1982.
13. Gusev A., Gerdt V. a. o. Symbolic-numerical algorithm for solving the time-dependent Shrodinger equation by split-operator method // Proc of the eigth Workshop on Computer Algebra in Scientific Computing, (September 12-16 2005, Kalamata, Greece). — 2005.
14. Сравнение алгоритмов для нормализации и квантования полиномиальных гамильтонианов / А. А. Гусев, Н. Чеканов, В. А. Ростовцев и др // Программирование. — № 2. — 2004. — С. 27-36.
15. LINA01: A REDUCE program for the normalization of polynomial Hamiltoni-ans / Y. A. Ukolov, N. A. Chekanov, A. A. Gusev et al // Computer Physics Communications. — Vol. 166. — 2005. — Pp. 66-80.
16. Puzynin I. V., Selin A. V., Vinitsky S. I. A high-order accuracy method for numerical solving of the time-dependent Schrodinger equation // Computer Physics Communications. — Vol. 123. — 1999. — Pp. 1-6.
17. Puzynin I. V., Selin A. V., Vinitsky S. I. Magnus-factorized method for numerical solving the time-dependent Schrodinger equation // Computer Physics Communications. — Vol. 126. — 2000. — Pp. 158-161.
18. Селин А. В. Метод приближенного решения линейного эволюционного уравнения в гильбертовом пространстве // ЖФМ и МФ. — № 42. — 2002. — С. 937-
19. Finite-element solution of the coupled-channel Schrodinger equation using highorder accuracy approximations / A. G. Abrashkevich, D. G. Abrashkevich, M. S. Kaschiev, I. V. Puzynin // Computer Physics Communications. — Vol. 85. — 1995. — Pp. 40-65.
20. Abrashkevich A. G., Kaschiev M. S., Vinitsky S. I. A New Method for Soling 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.
21. Wilcox R. M. Exponential operators and parameter differentiation in quantum physics // J. Math. Phys. — Vol. 8. — 1967. — Pp. 962-982.
22. Magnus W. On the exponential solution of differential equations for a linear operator // Commun. Pure Appl. Math. — Vol. 7. — 1954. — Pp. 649-673.
1963.
949.
e—
e
—ф 0
Вестник РУДН, Серия Математика. Информатика. Физика. № 1. 2008. с. 43-53 53
23. Crank J., Nicholson P. A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type // Proc Cambridge Philos. Soc. — Vol. 43. — 1947. — Pp. 50-67.
24. Березин И. С., Жидков Н. П. Методы вычислений. т.1. — М.: Физматлит, 1959.
25. http://www.nag.co.uk/numeric/numerical\_libraries.asp.
26. Калиткин Н. Н. Численные методы. — М.: Наука, 1978.
27. Физическая энциклопедия. т.5. — М.: Большая Россйская энциклопедия, 1998.
28. Yukalov V. I., Marzlin K.-P., Yukalova E. P. Resonant generation of topological modes in trapped Bose-Einstein gases // Phys. Rev. A. — Vol. 69. — 2004. — Pp. 023620-1-023620-16.
29. Olver P. J. Application of Lie groups to differential equations. — New-York, Tokyo: Springer-Verlag, 1986.
30. Обобщенный непрерывный аналог метода Ньютона для численного исследования некоторых нелинейных кванто-полевых моделей / И. В. Пузынин, И. В. Амирханов, Е. В. Земляная и др // Физика элементарных частиц и атомного ядра. — Т. 30. — 1999. — С. 210-268.
UDC 519.633.2, 517.958: 530.145.6
Multi-Layer Schemes for Solving the Time-Dependent Shrödinger Equation
O. Chuluunbaatar
Joint Institute for Nuclear Research 6, Joliot-Curie str., Dubna, Moscow Region, 141980, Russia
The algorithm based on unitary evolution operator decomposition to generate in a MAPLE and REDUCE packages multi-layer implicit schemes for numerical solving the time-dependent Shroedinger equation is presented. The optimal methods for construction of additional gauge transformations to extract symmetric operators needed for generation economical algebraic evolution schemes with respect to spatial variables by the finite element method are studied. The efficiency of the developed computational schemes till sixth order with respect to the time step and till seven order with respect to the spatial step, are demonstrated on the integrable model of oscillator in a time-dependent external field.
e—
e
—e e
