Диссертация (1150857), страница 5
Текст из файла (страница 5)
Приэтом интегрирование осуществляется только по области перекрывания сплайнов, которая довольно мала вследствие сильной локализации базисных сплайнов Эрмита.— 35 —Отметим, что разложение, в соответствии с формулами (1.60) и (1.55),каждой релятивистской компоненты биспинора по одинаковому наборуфункций (при этом конкретный вид этих функций не важен) приводитк наличию в спектре уравнения (1.58) нефизических (шпуриозных) состояний [59–61]. Однако, поcкольку вероятность перехода между физическим ишпуриозными состояниями достаточно мала, их наличие не оказывает влияние на расчёт процесса перезарядки.Рассмотрим нестационарный случай.
Запишем разложение (1.51) зависящей от времени волновой функции (, , , ) в следующим виде:(, , , ) =∑︁(1.66) () (, , ),=1где (,,) (, , ) ≡ ,(, , ) и единый индекс ≡ (, , ) определен-ным образом нумерует все базисные функции из набора (1.52). Подставивразложение (1.66) в нестационарное уравнение Дирака (1.25), получим систему уравнений, которую в матричной форме можно записать какS= H() ().(1.67)Здесь S – матрица перекрывания, H – матрица гамильтонина:′(1.68)S′ = ⟨ |′ ⟩ = ⟨,|′ ,′ ⟩,′H′ () = ⟨ |rot ()|′ ⟩ = ⟨,|rot ()|′ ,′ ⟩.(1.69)Более явные выражения для данных матричных элементов можно получить, если выполнить интегрирование по переменной , используя определение (1.52) базисных функций ,(, , ).
Для элементов матрицы пере-крывания получим, что′⟨,|′ ,′ ⟩ = B′ ′ ′ .(1.70)— 36 —Используя формулы (1.26) и (1.27), гамильтониан rot () можно записать как(1.71)rot () = () − () (Σ + ) ,где () – скорость вращения межъядерной оси, операторы (), Σ и всистеме координат (, , ) имеют вид (1.16), (1.28) и (1.30), соответственно.′Для того, чтобы получить ⟨,|rot ()|′ ,′ ⟩, найдем матричные элемен-ты для операторов в правой части выражения (1.71). Из того факта, чтогамильтониан () инвариантен относительно поворотов вокруг оси , и изформул (1.16) и (1.52) следует, что′⟨,|()|′ ,′ ⟩ = (H ()) ′ ′ , = 0 + , ′ = ′ 0 + ′ . (1.72)(, , ) и явный видПодставив формулу (1.52) для базисных функций ,оператора Σ в выражение для соответствующего матричного элемента, получим⎧⎨ ,′ +1 ,′ −1 , = 1, 31′⟨, |Σ |′ ,′ ⟩ = B′.⎩ ′ ′ , = 2, 42, −1 , +1(1.73)Аналогично, с помощью формулы (1.30) для оператора , найдем, что′⟨,| |′ ,′ ⟩(︂[︂(︂]︂′ )︂(−1)′Y′ − Z′ += ′ , ′ +1 X′ − +2(1.74)[︂(︂]︂)︂′ )︂(−1)′, ′ −1 −X′ − +Y′ + Z′ ,2где∫︁ X′ =′ ,∫︁Y′ = ′ ,∫︁Z′ = 2 (1.75)′ .Матричные элементы X′ , Y′ и Z′ рассчитываются численно, аналогичноэлементам (1.65)— 37 —Матрицы S и H() в рассматриваемом базисном наборе являются вещественными, что ускоряет вычисления.
Все матричные элементы, которые требуют численного интегрирования (см. выражения (1.65) и (1.75)), кроме V′ ,не зависят от времени и могут быть рассчитаны только один раз. Матрица Vзависит от времени через потенциал движущихся ядер, и её необходимо пересчитывать на каждом шаге по времени, поэтому в данном случае скоростьчисленного интегрирования является критическим параметром. Качественноускорить расчёт матрицы потенциала удается благодаря её разрежённости инебольшой области перекрывания БСЭ.1.3.4Решение стационарного уравнения ДиракаКак было показано в параграфе 1.3.3, стационарное уравнение Дирака в конечном базисе с помощью вариационного принципа сводится к обобщённойзадаче на собственные значения (1.58).
Используя стандартные численныепроцедуры, можно диагонализовать матрицу гамильтониана и получить всееё собственные вектора и собственные значения . Однако, полная диагонализация матрицы является очень ресурсоёмкой операцией и для достаточно большого базиса (∼ 104 ) требует слишком много времени и слишком большой объём оперативной памяти компьютера. При этом для расчёта процессаперезарядки необходимо получить лишь один собственный вектор, соответствующий начальному состоянию мишени.
Поэтому в представленной работеиспользовался описанный ниже подход, который позволяет быстро получитьвектор состояния с минимальной по модулю энергией (данным свойствомобладает основное состояние иона). Ранее подобный метод применялся длярешения стационарного уравнения Дирака в работе [20].Вначале найдём наибольшее по модулю собственное значение матрицы гамильтониана max и соответствующий ему собственный вектор max .
Для этой— 38 —цели применим итерационный степенной метод (см., например, [62]):+1S−1H ,= 0−1‖S0 H ‖(1.76)+1 = †+1 S−10 H +1 .Здесь ≥ 0 – номер итерации, вектор с увеличением будет сходитьсяк max , а будет стремиться к max . Начальный вектор 0 может быть выбран произвольным образом. Рассмотрим следующее матричное уравнениена собственные значения:[︁(︀ −1)︀2 ]︁2max − S0 H = 2 .(1.77)Наибольшее собственное число 2max и соответствующий ему собственный вектор в данном случае также могут быть получены с помощью степенного метода аналогично (1.76).
При этом будет выполняться условие2max = 2max − 2min ,(1.78)которое позволяет найти минимальное по модулю собственное значение minматрицы гамильтониана H , соответствующий ему собственный вектор будетсовпадать с таковым для 2max .Матрица H является разреженной, поэтому требует незначительный объём оперативной памяти для своего хранения, и умножение H на векторпредставляет собой достаточно быструю операцию, так как в неё вовлекаютсятолько ненулевые матричные элементы. Однако, матрица S−10 уже являетсяплотной, что вновь приводит к необходимости выделять слишком большойобъём оперативной памяти.
Кроме того, умножение S−10 на вектор будет наиболее затратной процедурой в ходе выполнения итераций степенного методаи окажет критическое влияние на скорость расчёта. Данную проблему можнорешить, используя тот факт, что матрица перекрывания S0 при соответствующей нумерации базисных функций представляет собой следующее прямое— 39 —произведение трёх матриц:S0 = E4 × S × S .(1.79)Здесь E4 – единичная матрица размерности 4 × 4, S и S – матрицы перекрывания одномерных БСЭ:(S ) ′max∫︁ () ′ (),=(S ) ′max∫︁= () ′ (),0(1.80)minгде () ≡ () и () ≡ () – одномерные кубические БСЭ, из которых строятся базисные функции (см. формулу (1.49)), индексы ≡ (, )и ≡ (, ) определённым образом нумеруют сплайны.
Из выражения (1.79) следует, что−1−1S−10 = E4 × S × S .(1.81)−1В результате надо хранить лишь очень небольшие матрицы S−1 и S , и нетникакой необходимости получать S0 −1 в явном виде. Таким образом, удаетсязначительно снизить требования к объёму оперативной памяти компьютера иодновременно существенно ускорить вычисления, так как разложение (1.81)даёт возможность оптимизировать умножение матрицы на вектор.1.3.5Решение нестационарного уравнения ДиракаМетод решения нестационарного уравнения Дирака или Шрёдингера с однойстороны должен обладать достаточным быстродействием, с другой стороныон должен гарантировать устойчивость и точность расчётов.
Более устойчивыми являются методы сохраняющие норму волновой функции, среди которых особо широко применяются метод Кранка-Николсон [63], метод сплитоператора [64, 65] и метод Арнольди-Ланкцоша [66]. Несмотря на низкую— 40 —стабильность, подходы, не гарантирующие сохранение нормы, также иногдаиспользуются [17, 20]. Целесообразность применения того или иного метода зависит от условий конкретной задачи.
Ниже описаны подходы, которыеиспользовались в данной работе.Формальное решение нестационарного уравнения Дирака в конечном базисе (1.67) может быть записано следующим образом (см., например, [67]):⎧⎫⎨ ∫︁out⎬−1(out ) = U (out , in ) (in ) = exp − S H( ) (in ).⎩⎭(1.82)inЗдесь in и out – начальный и конечный моменты времени, соответственно,U(out , in ) – матрица эволюции (пропагатор), и последнее равенство означает,что∫︁out∫︁out∫︁out∞∑︁(−)1 2 .
. . ×U (out , in ) =!=0in{︀inin(1.83)}︀S−1 H(1 )S−1 H(2 ) . . . S−1 H( ) ,где – оператор временного упорядочивания. Методы численного решениянестационарного уравнения (1.67) основаны на разбиении временного интервала [in , out ] на подинтервалы Δ , при этом матрица временной эволюциидля полного интервала принимает вид [67]U (out , in ) = −1∏︁U ( + Δ , ) ,(1.84)=0где 0 = in , +1 = + Δ и = out . Последовательно действуя пропагаторами U ( + Δ ) на начальный вектор (in ), можно получить конечныйвектор (out ).Далее, для упрощения выражений, по умолчанию будет предполагаться,что шаг по времени равномерный (Δ = Δ). Однако, необходимо отметить,что использование неравномерного разбиения интервала [in , out ] также возможно, и, более того, в определённых случаях путём оптимального выбора— 41 —такого разбиения можно улучшить баланс между быстродействием и точностью алгоритма.Пропагатор для интервала Δ можно записать как⎫⎧+Δ∫︁⎬⎨−1 S H( ) =U( + Δ, ) = exp −⎭⎩⎧⎫+Δ∫︁⎨⎬−1exp − S H( ) + (Δ3 ).⎩⎭(1.85)Последнее равенство в данном выражении следует из разложения Магнуса [68, 69].
Для достаточно короткого промежутка Δ можно считать, чтоH() не изменяется существенно в течение Δ, и применить следующую формулу:⎫⎧+Δ∫︁⎬⎨{︀}︀−1 S H( ) = exp − Δ S−1 H( + Δ/2) + (Δ3 ). (1.86)exp −⎭⎩Матрицу{︀}︀Ũ( + Δ, ) = exp − Δ S−1 H( + Δ/2)(1.87)будем использовать в качестве приближенного пропагатора:( + Δ) = Ũ(, + Δ) () + (Δ3 ).(1.88)Получить вектор ( + Δ) можно, диагонализовав матрицу H( + Δ/2).В её собственном базисе Ũ(, + Δ) также имеет диагональный вид.
Такойметод сохраняет норму волновой функции и является наиболее точным врамках приближения (1.88) для заданного шага Δ. Однако, диагонализацияматрицы является чересчур затратной процедурой, и на практике такой подход применим только при очень небольшом размере базиса (порядка сотнибазисных функций).— 42 —В качестве альтернативы можно воспользоваться разложением оператора Ũ по степеням Δ:]︂[︂)︀21 (︀ −1−1S HΔ + . . . (), ( + Δ) ≃ 1 − S HΔ −2(1.89)где H ≡ H(+Δ/2). Оборвав данное разложение на некотором члене, можнополучить (+Δ). Численный расчёт при этом сведётся к определенному количеству матрично-векторных умножений.
Умножение матрицы H на векторявляется относительно быстрой операцией вследствие разреженности H. Базисный набор, который используется в настоящей работе, даёт возможностьдля умножения вектора на матрицу S−1 применить следующее разложение:S−1 = E × S−10 .(1.90)Здесь E – единичная матрица размерности × , где – число каналов с различными проекциями углового момента на межъядерную ось, матрица S0 определяется выражением (1.63). C помощью подстановки разложения (1.81) для S−10 , формулу (1.90) можно записать в виде−1S−1 = E × E4 × S−1 × S .(1.91)Разложение (1.91) позволяет существенно ускорить умножение вектора наматрицу S−1 . Таким образом, алгоритм (1.89) характеризуется высокой скоростью работы на каждом временном интервале. Однако, серьезным недостатком такого метода является то, что он не сохраняет норму волновойфункции, приводя к неустойчивости расчётов.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
