Многослойные схемы для численного решения нестационарного уравнения Шрёдингера методом конечных элементов Текст научной статьи по специальности «Физика»

УДК 519.633.2, 517.958, 530.145.6

Многослойные схемы для численного решения нестационарного уравнения Шрёдингера методом

конечных элементов

О. Чулуунбаатар

Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6., г. Дубна, Московской обл., Россия, 141980

Построены симметричные операторно-разностные многослойные неявные схемы для решения нестационарного уравнения Шрёдингера до шестого порядка точности по шагу временной переменной на основе факторизации унитарного оператора эволюции с помощью явного разложения Магнуса. Выведены редуцированные схемы решения задачи Коши для системы связанных нестационарных уравнений Шрёдингера по гиперрадиальной переменной, используя разложение Канторовича волнового пакета по набору подходящих параметрических базисных угловых функций. Сформулированы неявные алгебраические схемы численного решения этой задачи с симметричными операторами, используя дискретизацию компонент волнового пакета по гиперрадиальной переменной методом конечных элементов высокого порядка точности. Эффективность и устойчивость построенных схем демонстрируется численным анализом точно-решаемых моделей одномерного осциллятора с частотой, зависящей от времени и двумерного осциллятора в переменном внешнем поле, используя стандартный базис угловых функций.

Ключевые слова: нестационарное уравнение Шрёдингера, задача Коши, факторизация оператора эволюции, операторно-разностные многослойные схемы, метод Канторовича, метод конечных элементов.

1. Введение

В работе [1] построены операторно-разностные многослойные (ОРМ) неявные схемы для решения нестационарного уравнения Шрёдингера (НУШ) на основе факторизации унитарного оператора эволюции через неявное разложение Магнуса [2], а также были сформулированы вопросы построения ОРМ схем с симметричными операторами, используя частичное расщепление оператора эволюции с помощью калибровочного преобразования, зависящего от функций в частном случае дипольного приближения в форме длины. Этот подход приводит к громоздким вычислениям коэффициентов разложения экспоненциального оператора эволюции с увеличением требуемого числа слоев ОРМ схем.

В данной работе развит более эффективный подход построения ОРМ схем с симметричными операторами на основе явного разложения Магнуса, т.е. разложения в ряд Тейлора логарифма оператора эволюции по шагу равномерной сетки временной переменной [2,3]. В этом подходе коэффициенты разложения оператора эволюции вычисляются в явном виде из системы рекуррентных соотношений. Используя частичное расщепление оператора эволюции с помощью калибровочного преобразования зависящего от операторов, построены ОРМ схемы с симметричными операторами до четвёртого порядка точности по шагу временной переменной для гамильтониана общего вида, а также шестого порядка точности при условии равенства нулю коммутатора первой и второй частных производных гамильтониана по временной переменной. Это обстоятельство позволяет использовать данный подход для решения более широкого класса эволюционных задач.

Статья поступила в редакцию 15 марта 2008 г.

Работа выполнена в рамках темы ОИЯИ «Математическая поддержка экспериментальных и теоретических исследований, проводимых ОИЯИ 09-6-1060-2005/2010» и при поддержке гранта РФФИ «Математическое моделирование динамики лёгких атомов и молекул под действием быстрых частиц, лазерных импульсов и магнитных полей 08-01-00604-а».

Содержание работы следующее. В разделе 2 представлен алгоритм вычисления оператора эволюции для решения задачи Коши НУШ и построена симметричная ОРМ схема решения задачи по временной переменной, необходимая для эффективного применения метода конечных элементов по пространственной переменной. В разделе 3 выведены редуцированные схемы решения задачи Коши для системы связанных НУШ по гиперрадиальной переменной, используя разложения Канторовича волнового пакета по набору подходящих параметрических базисных угловых функций. Сформулированы неявные алгебраические схемы численного решения этой задачи с симметричными операторами, используя дискретизацию компонент волнового пакета по гиперрадиальной переменной методом конечных элементов высокого порядка точности [4,5]. В разделе 4 стабильность и эффективность построенных алгебраических схем показана численным анализом точно решаемых моделей одномерного осциллятора с частотой, зависящей от времени, и двумерного осциллятора в переменном внешнем поле, используя стандартный базис угловых функций. В заключении обсуждаются основные результаты и дальнейшие перспективы применения построенных схем.

2. Общие формулировки и ОРМ схема

Рассмотрим задачу Коши для НУШ на интервале времени £ £ [¿о,Т] для начального состояния Ф(£о):

^^ = Н (¿)Ф(£), Ф(£о) = Фо, (1)

где Н(¿) — линейный самосопряжённый оператор. Потребуем непрерывность решения Ф(г,£) £ <8> о, Т]). Для решения задачи используется равномерная сетка по переменной £

Пт[¿о, Т] = {¿о, = 4 + т, (к = 0,1,..., К — 1), ЬК = Т} (2)

с шагом т. Решение Ф(£) в момент времени £ связано с начальным условием Фо с помощью оператора эволюции и(¿, ¿о)

(М°) = Н (¿)и (Мо), и (М) = 1. (3)

Решение Ф(£к+1) на к+1-м слое определяется через решение Ф(£к) на предыдущем к-м слое с помощью оператора эволюции и (¿к+1, ¿к) [6-9]

Ф(4+1) = и (¿к+1,£к)Ф(£к), (4)

и(¿к+1,£к) = ехр (-1тАкМ}(£к+1)) + О(т2М+1). (5)

Предположим, что оператор Н(¿) и его производные до порядка 2М — 1 являются гладкими функциями временной переменной ¿. Тогда, используя разложение Тейлора оператора Н(¿) в окрестности точки ¿с = ¿к + т/2

2М-1 (, _ , )*

Н(¿) = д*Н(¿с) + О (т2М) , (6)

*=о

разложим оператор АкМ) (¿) в ряд по т при £ = ¿к+1

2М-1

АкМ)(£к+1)= 5 т*Аш(£с),

3 = о

с неизвестными коэффициентами Заметим, что оператор и(^,^+1) яв-

ляется обратным к оператору и), из чего следует А^]^) = А^](^+1).

Так как выражение для а£М+1 ) получается из (7) формальной заменой то ряд (7) содержит только чётные степени т.

2.1. Ряд Тейлора для логарифма оператора эволюции

Рассмотрим уравнение (3) с эволюционным оператором и(¿, 8)

ди^ = -Н (*)и (¿,8), (8)

где 8 = ¿с — т/2 и £ = ¿с + т/2. Используя операторное тождество

= —и (М)и (М), и (8,£) = и-1 (¿,8), (9)

получаем следующее соотношение для оператора и(8, ¿)

ди (8,£)

т т,

«и (8,£)Н (¿). (10)

Если заменим £ ^ 8 в уравнении (10), то приходим к соответствующему операторному уравнению для и(¿, 8)

ди|;^ = ,и (¿,8)н (8). (11)

Отсюда, используя уравнения (8), (11) и равенство

2 ди (¿,8) = ди (¿, 8) ди (¿, 8)

дт д8 , ()

получаем выражения для производной оператора и(¿,8) по параметру т при фиксированном ¿с

2ди(^ 8) = — гЯ(¿)и(¿, 8) — ги(¿, 8)Н(8). (13)

дт

Используя известную формулу [10] для вычисления производной экспоненциального оператора и(¿, 8) = ехр(Р) при Р = — гтА^] (¿)

1

ди(¿,8) _ > д^

дт

0

Г дР

ёж ехр(жР)— ехр(—жР )и (¿,8), (14)

и формулу Кампбелла-Хаусдорфа [3]

ехр(А)В ехр(—А) = Е 1(ойА)3' В, (15)

5=0 - '

перепишем уравнение (13) в следующем виде

1 дР _ 1

2 Е ИРУ ^ = — (¿) — 1 £ - (а^РН(8). (16)

Здесь линейный оператор (а^А) : £(Х) ^ £(Х) (С(Х) — пространство линейных

операторов) определён для операторов А, В £ С(Х) в виде (а^А)В = [А, В] =

AB — BA и обладает следующими свойствами: (adA)0B = B, (adA)jB = (adA)j 1 (adA)B. Используя разложения (6) и (7) получаем

(_9)qтq / \ q тi

— Е Нг- Е (ic) Е тnQqn(tc) = Е ^Тм dtH(^

q=0 \ m=0 ) n=0 i=0 (17)

( — 1)" П + 1 -Qqn(tc) = 2П+4 d"H (ic) — ^ А(„)&).

Тогда, неизвестные коэффициенты вычисляются в явном виде из системы рекуррентных соотношений

1

+ 1)А(^С) = 2*+1./ Н(^ + ЕЕ Е -, (18)

■>' д=1 п=о г1,...,г,=о 4'

(¿с) = МА^^))... МА^ ^с))^^).

При М = 1, 2, 3 операторы АкМ) = АкМ) + гА1кМ) = 4М)(^к+1) имеют вид

= Н, 4^ =0, 2 2 42) = 44к1) + 24 Н, 44к2) = 4^ + ^МН) Н,

т4 .... т4 .. т4 . (19)

АГ = Ak2) + — H — Т-(adH)2 H — -V(ad H)2H,

k k 1 non >70 n 4 ' олпч ' '

{(3) _ h — (adH)2 H — — ^

1920 720 ; 240

444

А3) = ^k2) — 4T80(ad H)H + ^(ad H) H +720(adH)3 H,

где H = H (ic), H = dtH (i)|t=tc ,...

2.2. ОРМ неявная схема

Для генерации операторно-разностной неявной схемы представим экспоненциальный оператор (4), (19) на каждом k-м слое в виде обобщённого [M/M] Паде разложения с той же точностью [6, 7]

exp (—irAf)) = П Tfk + O(r2MТ1),

С=М

Т + та[М ЧМ) \-1 + т«'М '4кМ) ^ (20)

^ = I7 + 2М J + 2М

где 7 — единичный оператор. Здесь коэффициенты а^М) (£ = 1,..., М, М ^ 1) являются корнями полиномиального уравнения 1Т\( —М, —2М,2Мг/а) = 0, где 1^1 (а, 6, ж) — вырожденная гипергеометрическая функция, коэффициенты ) комплексно сопряжены с а^М). Коэффициенты а^М) имеют следующие свойства [6]: ) < 0 и 0,6 < |а[М)| < где ^ ~ 0, 28 корень трансцендентного уравнения ^ехр(^ + 1) = 1. Заметим, что условие

т< 2М^||4кМ)(^)||-1, (21)

гарантирует сходимость разложения (20) для ограниченного оператора А^

Используя разложение (20), произведём факторизацию эволюционного оператора, т.е. представим (5) в виде системы М уравнений относительно набора М — 1

вспомогательных функций

^0 = ),

та(м)А(м)\ , ( та(м)А(м)\ , „

I + ^^^ ) ^ = и + ) , С =1,...,М, (22)

Ф(^+1) =

решение которой обеспечивает переход от ) на к-ом слое к Ф(^+1) на к + 1-ом слое. Заметим, что в данном случае сохраняется унитарность приближенного оператора перехода, поскольку оператор ) самосопряжённый. Поскольку ) = 0, вспомогательные операторы Т<^ изометрические, следовательно норма всех ||^/М || равны

'•0" = II = ••• = ИЛ (23)

2.3. Симметричная ОРМ неявная схема

Для получение схем с частичным расщеплением оператора эволюции, т.е. с выделением симметричной части )(^) оператора )(^), применяется калибровочное преобразование = ехр (гБ^

аА£М) (¿) = ехр (4м) (*)) А<м) (¿) ехр (—4м) (*)) . (24)

Оператор ) (¿) ищется в виде ряда по степеням т при ^ =

2М-1

^М)(4+1)= Е т5%)(*с), (25)

5=0

с коэффициентами которые определяются из условия

А1М)= ехр (4м^+1)) А<м) ехр (—^+1)) = 0(т2м). (26)

Поставляя разложение ) = ) (¿^+1) в уравнение (26) и приравнивая коэффициенты при одинаковых степенях , получаем систему рекуррентных уравнений для нахождения искомых коэффициентов = б^)^) с начальным условием £(0) = 0. Первые три оператора (25) и (24) имеют вид

2.

^ =0, £(2) = £(1) + ^ Н,

4 4

£(3) = ^^ + 480 Н +72о(а^Я)2 Я, при МЯ) Н = 0,

(27)

2..

А^ = а!£1) = Н, А12) = А12) = АА£1) + 2, Н,

4 . 4 .... 4 .. 4 .

А[3) = А[3) + — (ай Н)2Н = А!2) + — я' — ^(айН)2 Н —Г—(ай Н)2Н. ь ^ 288 1920 720 ^ ; 1440 ;

В результате, на каждом к-м слое сетки Пт[¿0,Т] (к = 0,1,. схему с частичным расщеплением оператора эволюции

, К — 1) получаем

I +

та

_(м )А/(м)

С

к

2М

= ехр (4м^ Ф^),

лС/м_ |7 + таГ)

^¡Г =

2М

Йс-1)/м, С = 1,

., М, (29)

Ф(4+1) = ехр ( —

?(м Р

^.

Эта унитарная М-слойная неявная схема сохраняет норму (23) разностного решения и является стабильной [11]. Каждое уравнение по отдельности из системы М уравнений (29) имеет порядок аппроксимации не выше второго, тогда как вся схема (29) имеет суммарную аппроксимацию 2М-го порядка О(т2м) по временному шагу сетки [11]. ( )

Далее представим экспоненциальный оператор ехр "*) на каждом к-м

слое в виде обобщённого [Ь/Ь] Паде разложения аналогично (20). Эта аппроксимация имеет порядок О (т4-+2), поэтому чтобы выполнялось условие 4Ь+2 ^ 2М, было выбрано Ь = [4^], где [ж] обозначает целую часть значения ж. В этом случае получаем модифицированную схему

'

_(м )А/(м)

I +

та

С

к

2М

10 = ),

10 = 11, (м )А/(м)

4п-1)/-,

п = 1,...,Ь,

= |1+

а

с

к

2М

Йс-1)/м,

С = 1,

М, (30)

10 = 11,

^ _(-)£(м) \ / а(-)£(М)-

I + аП-2Ьм^ = I + ап

4П-1)/Ь,

п = 1,...,Ь,

2Ь

фс^+1) = 11.

3. Применение симметричной ОРМ неявной

схемы для НУШ

Рассмотрим й-мерное НУШ с гамильтонианом Н(¿) = Н(г,^) д Ф(г,<)

д*

= Н (г, ¿)Ф(г, ¿), Ф(г, ¿0) = Ф0(г),

1

(31)

Н (г, ¿) = Н0(г) + f (г, ¿), Н0(г) = — -У2 + и (г).

Потребуем непрерывность волновой функции Ф(г,<) е Wl(Rd ® [¿0,Т]). Также предполагаем, что эти решения описывают модели атома во внешнем поле f (г,^) с f(г, ¿0) = 0, и функция f (г, ¿) имеет непрерывные частные производные по временной переменной до порядка 2М 1 и по пространственной переменной до порядка 2М — 2. Волновая функция Ф(г,^) удовлетворяет условию нормировки

||ФСМ)||2 = I |Ф(г,*)|2ёг = 1, г е [¿0,Т].

(32)

Первые три оператора 4М) и 5*кМ) имеют вид

4^ = Н, ^к1) = 0,

42) = 4^ + с(2), ^к2) = ^к1) + 2(2), (33)

4 4

43) = 42) + ^(3) — 7^V, (У2 /) V,, ^3) = *(2) + 2(3) + V, (У2 /) V,,

где функции С(2), С(3), Z(2) и Z(3) определяются в терминах / = /(г,^с), /=

д/ (г,^ ,..., и = и (г):

2 2 с(2) = — /', 2(2) = — /,

^ 24 ^ 12 ],

4 4 2 4 4

с(3) = 1920 7' +1440 - 720 ') (^<и + />) — 2880 (V4 '

4 4 4

2(3) = 480 ' ') (^<и + /0 + 5880 ')•

(34)

Далее используется схема только при М ^ 3, поскольку при М ^ 4 рассматриваемая схема содержит оператор набла с третьего порядка.

3.1. Разложение Канторовича

В методе Канторовича [12] парциальная волновая функция Ф(г,^) ищется в виде разложения по однопараметрическому набору гиперповерхностных функций {В* (П; г)}*=а1Х:

*тах

Ф(г,*) = 5 В*(П;г)х*М). (35)

*=1

В разложении (35), вектор-функция х(г, ¿) = (%1(г, ¿),..., Х*тах(г, ¿))Т является неизвестной, а поверхностная вектор-функция В(П; г) = (В1 (П;г),...,В*тах(П;г))Т зависит от набора угловых переменных П для каждого значения гиперрадиуса г, который рассматривается здесь как параметр. Базисные функции В*(П;г) £ ~ Ь2($^-1(П)) определяются как решения параметрической задачи на собственные значения

(-Л2 + и(г)) В*(П;г) = Я,-(г)В*(П;г), (36)

где Л2 — самосопряжённый оператор обобщённого углового момента. Собственные функции этой задачи удовлетворяют краевым условиям по угловым переменным П, аналогичным тем, которые наложены на волновую функцию Ф(г,^) при каждом фиксированном значении параметра г £ Н+ и условию нормировки

^ВДП; г) В*(П;г)^ = ^ ВДП; г)В*(П;г)ёП = %, (37)

где — символ Кронекера, а черта обозначает комплексное сопряжение.

После вариации функционала Рэлея-Ритца с использованием разложения (35) уравнение (31) сводится к системе из .тах обыкновенных дифференциальных уравнений второго порядка относительно неизвестной функции х(г, ¿):

1! дх(М = Н(г, ¿)х(г, ¿), х(г, ¿о) = Хо(г),

Н(м) = - ^1 ¿г^-1 £ + V(г) + 2 (1)М) + д(г) ^ + ^

Здесь I, У(г), Ъ(1)(г, ¿) и Р(г) матрицы размерностью ^тах х ^тах, элементы которых определяются соотношениями

^ (г)= Е«(г) + Е(г) . + 1/ дВ«(П; г) ^ (г)= 2 + 2\ дг

■5г(,1)(г,¿) = <В«(П;г) f(г,¿) В^-(П;г)),

дВ,- (П; г)

дг

^ (г) = — В«(П; г)

дВ5- (П; г)

дг

Граничные условия для функции х(г, ¿) имеют вид

х(0, ¿) = 0, если 11ш г^-115(г, ¿)| = то,

Г >0

(39)

11ш г^ 1 — р(г)^х(г, ¿) = 0, если ш1п 11ш г^ (г, ¿)| < то, (40)

11ш х(г, ¿) = 0.

Функция х(г, ¿) удовлетворяет условию нормировки в соответствии с формулой (32)

У (Х(г, ¿))Т х(г, ¿)г^-1ёг = 1.

(41)

В этом случае получаем матричную операторно-разностную схему для неизвестных векторов х(г,¿), аналогичную схеме (30)

^ _(-)/ (м )■

I _ ап °к

2Ь

~п/-X к

х к = хМк^

/ а(-)5 (м )'

/1 _ ап » к

2Ь

01 X к = X к,

X к

(Т1)/-, П =1,

Ь,

(м )Л (м )\ 1 + тас Лк | с/м

2М 1 ^

I +

та[м км Г 2М

X Г1)/м, С = 1,...,М, (42)

( _(-)й (м)' 1+^

~п/-X к

х к = х: к,

^ а(-)5 (м )'

I + ап ° к

~(п-1)/- 1

Xк' )7 , V = 1,

2Ь

х(г,гk+l) = х: к.

Здесь при М = 1, 2, 3 операторы Л¿м) и 5¿м) имеют вид

.,Ь,

(1)

Н(г, ¿с), 5к1) = 0,

Лк2) = Лк1) + с(2), 5к2) = 5к1) + ъ(2),

,(з)

Л к2) + с (3)+ ск3),5

. (3)

(3)

(2)

+ ъ(3) — с

(3)

(43)

где (С(1) и Ъ(1), I = 1,2, матрицы размерностью ^тах х ^тах, элементы которых определяются соотношениями

с^Мо)=с?.

= п « -

В«(П; г) С(1) В,-(П; г) ^(гЛ) = = (в«(П; г)|- (1)|В5 (П; г))

п

п

к

к

п

п

Самосопряжённый оператор Ск3) имеет вид

С

(3)

1 д

720 \ 1 дг

д

г^-10(г) — + Л(г) - ((т(г)^- +

дг -1-

д 1 дг^-1С((г)'

дг г^ 1 дг

, (45)

где 13(г,¿с) = 1)(г), Л/"(г, ¿с) = У(г) и С((г, ¿с) = (((г) — матрицы размерностью ^тах х ^тах, элементы которых определяются соотношениями

Я,(г) = <Бг(П;г) V? л В(П;г)

% (г)

дД(П; г)

дг

V? /

дБ,-(П; г) \ (Б«(П; г)|д(г)|Б,(П; г))

^ / + о

дг

<9, (г) = -( Бг(П; г)

V? /

дБ, (П; г)

(46)

дг

£(г) = г2 [Е(г) + Е(г) - 2и(г)] (V? Л + 1 (л2, (V? /

4

к

2

2

г

3.2. Приближения высокого порядка в МКЭ

Для численного решения задачи (38) на сетке Пт[¿о,Т] граничные условия (40) и условие нормировки (41) по пространственной переменной г на бесконечном интервале заменяются соответствующими условиями на конечном интервале Пг [гт1п, гтах]. Далее, на каждом к-ом шаге сетки Пт [¿о, Т] рассмотрим дискретное представление решений х(г, ¿к) задачи (38) с помощью метода конечных элементов на сетке = {гтт = го, г, = г,— + Л, , гп — гтах } в виде конечной суммы локальных функций ЛР(г)

ХМ(г, ¿к) = Е ) Л? (г), ^ — 1,... ,^тах, (47)

?=о

где Х^(^к) являются узловыми значениями неизвестной функции хДг, ¿к), относительно которых численно решается исходная задача. Локальные функции Л?» являются кусочно-непрерывными полиномами порядка р [13]. Коэффициенты ) формально связаны со значениями вектора задачи хДг, ¿к) в Лагран-жевых узловых точках [13]

Л •

гр,г = г,—1 +—- г, Л, — г, - г,—1, г — 0, ...,р, (48)

р

соотношениями

Х^к) = Х^р^к), I = г + р(? - 1). (49)

Подставляя разложение (47) в операторно-разностную схему (42), умножая слева на ^р(г) и интегрируя на интервале Пг[гт1п,гтах] рассматриваемая схема (42) приводится к системе алгебраических уравнений для вектора Хк =

Г , „ Ч пр 4 5тах

¿к)о| ^ при заданном М

х к = Xk,

(вр — а^8^ Xк/- = (вр — р) хкп-1)/-, п = 1,...,Ь,

х к = х к, _(м) \ / (м)

вр + ^аМ-Лк I хk/м = (вр + ^Лр) Xг1)/м, С = 1,...,М, (50)

х к = х к,

вр+а^вк! хк/-=^вр+а^в ^ хг1)/-, п=1,...,ь,

хк+1 = х к.

Симметричные матрицы Лр, вр и Б к имеют ленточную структуру, кроме того матрица вр положительно определена. Они имеют вид

Лк = Е ар,к, вр = £ ьр, Бк = Е 5. (51)

5 = 1 5 = 1 5 = 1

где локальные матрицы ар к, Ьр и 8р к вычисляются по формулам, аналогичным представленным в работе [13]. Эти интегралы вычисляются с помощью Гауссов-ской квадратурной формулы порядка р+1. Как известно [4], теоретическая оценка

по норме || ■ || на конечно-элементной сетке П^ между точным х(г,¿к) и приближенным хк решениями имеет порядок

11х(г, ¿к) — хк II = О(^р+1), (52)

где Н — максимальный шаг сетки, так как можно выбрать т < Н и р ^ 2М согласно оценке (21).

Для анализа сходимости предложенных схем на последовательности трёх вдвое сгущающихся временных сеток вычислялись абсолютные погрешности

тах

Ег^Н Е У X(г, ¿) — XV3' (г, ¿)|2И-Чг, з = 1, 2, 3, (53)

^=1 0

и коэффициент Рунге [14]

= 1082

Ег^, 1) — Ег^, 2)

Ег^, 2) — Ег^, 3)

(54)

где XV3 (г, ¿) — численное решение с шагом по временной переменной т5- = т/25-1. В качестве xV(г, ¿) используются либо точное решение, либо численное решение с шагом т4 = т/8. Отсюда получаем численные оценки для порядка сходимости предложенной схемы (50).

4. Эталонные модели

4.1. Одномерный осциллятор

НУШ для одномерного осциллятора с переменной частотой на конечном интервале £ £ [0, Т] имеет вид

1 _о^ (£)ж2

2 дж2 + 2

£), Мж) = ехр (- 1(ж -^2)2)

(55)

где ш2(£) =4 — 3ехр(—£) [6]. Точное решение задачи Коши для уравнения (55) имеет вид

^0X1 (ж, £) = ехр (-X(£)ж2 + 2У(£)ж — 2(£))

(56)

где функции X(£), У (£) и 2(£) удовлетворяют задаче Коши для системы уравнений

1 ¿х(£) = 2х2(£) , х(0) = 2, г Ау (£) = 2Х (£)У (£), У (0) = ^,

1 (£) = -X (£) + 2У 2(£), 2 (0) = 1.

(57)

Для аппроксимации решения по переменной ж использовались конеч-

ноэлементная сетка Пш[жт;п,жтах] = (жт;п = —10(100) жтах = 10}, где число в скобках означает число конечных элементов на интервале, и шаг по времени т = 0, 009765625. Между каждыми двумя узлами сетки применялась интерполяционная формула Лагранжа восьмого порядка (р = 8). Проведены численные эксперименты на сгущающихся сетках ПЖ[жтш, жтах], которые показали строгое соответствие коэффициента Рунге (54) теоретическим оценкам (52) при р = 8.

На рис. 1 показаны абсолютные погрешности Бг(£;] = 1, 2, 3, для схемы с М = 1, 2, 3. Возле каждой тройки кривых показано среднее значение вм (£) по всем значениям в(£к) на сетке Пт [0,10]. В результате подтверждены теоретические оценки порядка сходимости предложенных схем (50), так как теоретический порядок равен в(£) = вм (£) = 2М. Расчёты выполнялись с двойной точностью 16 значащих цифр для М = 1, 2. Для М = 3 использовалась четверная точность 32 значащих цифр, поскольку абсолютная погрешность Бг(£; 3) меньше чем 10-14 в окрестности £ = 0.

10"' 10* 10л 10" 10* 10-*

=5101 е. ю

Ш 10*10" 10"11 10" 10"13 10""

2М=2 р,=1.99

2М=4 р,=3.9Е

2М=6 =5.99

123456789 10

I

Рис. 1. Абсолютные погрешности Ev(t,j), j = 1, 2, 3 (штрих-пунктирная, пунктирная и сплошная кривые), соответственно для схем второго, четвёртого и шестого порядков с шагом по времени т = 0, 009765625.

4.2. Двухмерный осциллятор

НУШ двухмерного осциллятора (или заряженной частицы в постоянном однородном магнитном поле) во внешнем электрическом поле Е (¿) = (Е1 (¿), ¿2^)) на временном интервале е [0, Т] в дипольном приближении имеет вид [15]

д ¡t

д2

д2

дж2 + ду2

гш f д

+ Т Г ду1

У1

дж1

^2 \ +у ^ + у2) — ж^ф — у^^ Ф(Ж1,У1,*). (58)

Во вращающейся с частотой ш/2 системе координат ж1 = ж 008(^/2) + у 8т(^/2) и у1 = у 008(^/2) — ж 8т(^2), уравнение (58) принимает вид

^ == (-2 (12+£)+т <*2+y2)+-/i (t)+y/2 y, t)

(59)

где /i(t) = -Ei (t) cos (wí/2) + £2 (t) sin (wi/2), /2ft) = -Ei(t) sin (wt/2) -£2^) cos (wt/2). В полярных координатах (r, 0) уравнение (59) сводит к уравнению

^(r,0,t) _ 1 f 1 + — —^

r--+ 1

2 V r дг дг r2 д02 У

+

ш2г2

8

+ r(/1 (t) cos(0) + /2 (t) sin(0)) ф(г, 0,t). (60)

Разлагая решение по угловым функциям Bj (0) G L2(S 1(0)) (метод Галёркина)

jmax

B1(0)

0(r,0,í)= £ Bj(0)Xj (r,t),

j=1

1 = , B2j (0) = si^j) , B2j + 1(0) = CO^^ , j > 1,

(61)

\/2п' " ' ' ^ " ' ' ' ^ получаем систему обыкновенных дифференциальных уравнений (38) относительно неизвестных коэффициентов (х5- (г, ¿)}5-=г1х на интервале времени t е [0, Т]

г I

/1 д д \ = (-2r1 дrrдr + ^J x(r,¿),

lim r дx(r, t) = 0, lim ^(r, t) = 0,

r—Q дr r—

(62)

где У(г, ¿) — семидиагональная матрица размерностью з'тах х з'тах, элементы которой определяются соотношениями

ш2г2

1

"8 + 2Т2

г = j

Vij (r, t) =

—2min(i,j) = 1, max(i,j) = 2, r

—2 min(i,j) = 1 max(i,j) = 3,

r (63)

-/2(í), mod(min(i,j), 2) = 1 и max(i, j) = min(i, j) + 1,

r

- -/2(t), mod(min(i, j), 2) = 0 и max(i,j) = min(i, j) + 3,

r

2 max(i,j) = min(i,j) + 2, 0, иначе.

i

2

Функции X] (г, £) в момент £ = 0 выберем в виде

/ 1

X](Г 0) = Vйехр 4] (64)

что соответствует функции основного состояния свободного осциллятора (58)

(ж,у) = \/1П ехр (-1 и(х2+у2)) ■ (65)

В этом случае точное решение задачи Коши для уравнения (59) с начальным условием (65) имеет вид

фе*(®,У,*)= ехр (-Х1(£)ж2 +2У1(£)Ж - ^(£)) х V 2п

х У2Пехр (-Х2(£)у2 + 2ВДу - , (66)

где функции X](£), у(£) и (£), j = 1, 2, удовлетворяют задаче Коши для системы уравнений

2

1 -X](£) = 2Х|(£) - , X](0) = 4,

1 Ау, (£) = 2Х] (£)У] (£) + М, У] (0) = 0, (67)

1 ^ (£) = -X] (£) + 2У/(£), (0) = 0.

Заметим, что X](£) = и/4. Эта задача имеет аналитическое решение в случае внешнего поля Е](£) = й] 8ш(и]£), ^ = 1, 2, что даёт тестовый пример для исследования эффективности численных алгоритмов и сходимости по числу ^тах радиальных уравнений. Проекции точного решения уравнения (60) по базису угловых функций (61) определяются по формуле

2п

ХГМ) = / В (0)фе*(г,М№0. (68)

В расчётах были выбраны следующие параметры: и = 4п, и>1 = 3п, и2 = 5п, й1 =24 и й2 = 9. При этом плотность вероятности |ф(г, 0,£)|2 = |ф(г, + Т)|2 является периодической функцией с периодом Т = 2.

Для аппроксимации решения X] (г, £) по радиусу г использовалась конечно-

элементная сетка Пг [гт;п,гтах] = (гт;п = 0(120)1.5(60) гтах = 4} и шаг по времени т = 0,0125. Между каждыми двумя узлами сетки применялась интерполяционная формула Лагранжа восьмого порядка (р = 8). На рис. 2 показаны абсолютные погрешности Ег(£; j), для схемы с М =1, 2, 3 при ^тах = 20 и ^тах = 30. Возле каждой тройки кривых показано среднее значение вм (£) по всем значениям ) на сетке Пт [0, 2]. В результате подтверждены теоретические оценки в(£) = вм (£) = 2М порядка сходимости предложенных схем (50). Как видно из рис. 2, при ^ах = 20 на интервале времени £ £ [0, 2] обеспечивается абсолютная погрешность решения Ег(£; 3), меньшая чем 10-1, 10-3 и 10-5, соответственно для М =1, М = 2 и М = 3,и при увеличении ^ах уменьшается.

На рис. 3 показаны абсолютные погрешности фехъ(ж, у, £) - ф(х, у, £) при £ =1, 2 и £ = 2. Здесь ^ах = 20, Т4 = т/8 = 0,0015625 и М = 3.

На рис. 4 показаны абсолютные погрешности для ,?тах = 30.

В работе [16] по каналированию и сверхфокусировке пучков лёгких ядер в тонких плёнках было указано, что главное свойство гармонического осциллятора —

t t

(a) jmax = 20 (b) jmax = 30

Рис. 2. Абсолютные погрешности Бг(£, ,;), М = 1, 2, 3 (штрих-пунктирная, пунктирная и сплошная кривые), соответственно для схем второго, четвёртого и шестого порядков с шагом по времени т = 0, 0125

(a) (b)

Рис. 3. Абсолютные погрешности фехь(х,у,%) — Ф(Х,У,£) при £ = 1, 2 и £ = 2. Здесь

jmax = 20, т4 = т/8 = 0, 0015625 и М = 3.

(a) (b)

Рис. 4. Тоже самое, что на рис. 3 при jmax = 30.

это повторная фокусировка волнового пакета в координатной и импульсной плоскостях. Как известно, идеальная линза преобразует плоскость изображений в её Фурье-образ. Согласно этому первоначально широкий гауссовский волновой пакет преобразуется в чрезвычайно узкий волновой пакет в фокусе. Это «дыхание»

пучка внутри канала, предсказанное по квантовомеханической модели двумерного осциллятора без внешнего поля, было подтверждено компьютерными экспериментами [16]. Для иллюстрации на рис. 5 приведена динамика замкнутой кривой, на которой плотность вероятности волнового пакета |ф(ж,у,£)|2 в рассмотренной точно решаемой модели двумерного осциллятора во внешнем поле в данный момент времени £ в два раза меньше максимальной. Расчёты проводились при и = 4п, и1 = 3п, и2 = 5п и й1 = й2 =0 (а), и й1 = 24, й2 = 9 (б). Волновая функция в начальный момент времени ф0(х,у) = и/(20п) ехр(-и(ж2 + у2)/40). Чёрным цветом показана динамика волнового пакета основного состояния свободного осциллятора (65).

Рис. 5. Динамика замкнутой кривой, на которой плотность вероятности |^(х,у,£)|2 в данный момент времени £ в два раза меньше максимальной:

|ф(х,у,£)|2 = шахж,у |ф(х,у,£)|2.

5. Заключение и обсуждение результатов

В работе представлен новый подход решения задачи Коши для НУШ, в котором используется симметричные ОРМ неявные схемы до шестого порядка точности по шагу равномерной сетки временной переменной и их редукция по пространственной переменной с помощью разложения волнового пакета проекционными методами Галёркина и Канторовича, а также МКЭ. На модельных расчётах явно показано строгое соответствие теоретическим оценкам численных оценок сходимости сгенерированных схем второго, четвёртого и шестого порядков точности относительно шага равномерной сетки по временной переменной и соответствие восьмого порядка пространственной переменной. Разработанный подход позволяет учесть влияние ангармонических поправок на динамику волнового пакета в модели двумерного ангармонического осциллятора во внешнем поле, а также изучать роль ангармонических поправок в задаче динамики двумерных систем под влиянием внешнего магнитного поля и термостата [17]. Применение ОРМ схем к расчёту модели трёхмерного атома водорода в магнитном поле при воздействии последовательности сверх коротких импульсов, используя параметрический базис угловых сплюснутых сфероидальных функций для разложения волнового пакета методом Канторовича, будет рассмотрено в следующей работе.

Автор благодарит С.И. Виницкого, А.А. Гусева, М.С. Касчиева, И.В. Пузы-нина и Л.А. Севастьянова за сотрудничество и поддержку.

Литература

1. Чулуунбаатар О. Многослойные схемы для численного решения нестационарного уравнения Шредингера // Вестник РУДН: Серия Математика. Информатика. Физика. — № 1. — 2008. — С. 59-69.

2. Magnus W. On the Exponential Solution of Differential Equations for a Linear Operator // Commun. Pure Appl. Math. — Т. 7. — 1954.

3. Wilcox R. M. Exponential Operators and Parametr Differentiation in Quantum Physics // J. Math. Phys. — Vol. 8. — 1967. — Pp. 962-982.

4. Strang G., Fix G. J. An Analysis of the Finite Element Method. — New York: Prentice-Hall, Englewood Cliffs, 1973.

5. Bathe K. J. Finite Element Procedures in Engineering Analysis. — New York: Englewood Cliffs, Prentice Hall, 1982.

6. Puzynin I. V., Selin A. V., Vinitsky S. I. A High-Order Accuracy Method for Numerical Solving of the Time-Dependent Schrodinger Equation // Comput. Phys. Commun. — Vol. 123. — 1999. — Pp. 1-6.

7. Puzynin I. V., Selin A. V., Vinitsky S. I. Magnus-Factorized Method for Numerical Solving the Time-Dependent Schrodinger Equation // Comput. Phys. Commun. — Vol. 126. — 2000. — Pp. 158-161.

8. Селин А. В. Метод приближённого решения линейного эволюционного уравнения в гильбертовом пространстве // ЖВМ и МФ. — Т. 42. — 2002. — С. 937949.

9. Селин А. В. Метод аппроксимации эволюционных операторов с помощью экспоненциального представления и рациональных функций в гильбертовом пространстве. Дисс. канд.физ.-мат. наук. — Дубна, 2002.

10. Snider R. F. Perturbation Variation Methods for a Quantum Boltzmann Equation // J. Math. Phys. — Т. 5. — 1964. — С. 1580-1587.

11. Самарский А. А. Теория разностных схем. — М.: Наука, 1977.

12. Канторович Л. В., Крылов В. И. Приближённые методы высшего анализа. — М.: Гостехиздат, 1952.

13. 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. — Vol. 177. — 2007. — Pp. 649675.

14. Калиткин Н. Н. Численные методы. — М.: Наука, 1978.

15. Butkovskiy A. G., Samoilenko Y. I. Control of Quantum-Mechanical Processes and Systems. — Dordrecht Hardbound, Kluwer Academic Publishers, 1990.

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

17. Kalandarov S. A. et al. Influence of External Magnetic Field on Dynamics of Open Quantum Systems // Phys. Rev. E. — Vol. 75. — 2007. — Pp. 031115-1-16.

UDC 519.633.2, 517.958, 530.145.6

Multi-Layer Schemes for Solving Time-Dependent Schrodinger Equation by Finite-Element Method

O. Chuluunbaatar

Joint Institute for Nuclear Research Joliot-Curie str., 6 Dubna, Moscow Region, Russia, 141980

The symmetric implicit operator-difference multi-layer schemes for solving the time-dependent Schrodinger equation based on decomposition of the evolution operator via the explicit Magnus expansion up to the sixth order of accuracy with respect to the time step are presented. Reduced schemes for solving the Cauchy problem of a set of coupled time-dependent Schrodinger equations with respect to the hyperradial variable are devised by using the Kantorovich expansion of the wave packet over a set of appropriate parametric basis angular functions. The implicit algebraic schemes for numerical solving the problem with symmetric operators, using discretization of the component of the wave package by hyperradial variable by the high order finite-element method are formulated. The convergence and efficiency of the numerical schemes are demonstrated in numerical calculations of the exactly solvable models of one-dimensional oscillator with time-dependent frequency, two-dimensional oscillator in time-dependent external field by using the conventual angular basis.