CC BY
35
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
СРАВНЕНИЕ ДВУХ АЛГОРИТМОВ ВЫЧИСЛЕНИЙ В РАЗНЫХ ФАЗАХ МОДЕЛИ ЭЛЕМЕНТАРНЫХ ЧАСТИЦ
Н.В. ЗВЕРЕВ, доц. каф. физикиМГУЛ, канд. физ.-мат. наук
Одной из важных задач квантовой теории поля является создание и реализация эффективных алгоритмов для исследований моделей элементарных частиц. Практический интерес представляет векторная модель частиц, построенная на матричной группе преобразований U(1) для обеспечения необходимой инвариантности действия модели относительно преобразований волновых функций частиц. Для аналитического и численного исследований моделей частиц Вильсоном предложен математический метод моделирования - метод решетки [9]. Этот метод основан на аппроксимации непрерывного пространства-времени дискретной совокупностью точек с заменой непрерывных волновых функций дискретными величинами. При исследовании моделей на решетке пространства-времени вычисляют их корреляционные функции, используемые для вычисления физических характеристик элементарных частиц. Вычисления корреляционных функций выполняют, главным образом, методом континуального интеграла по всем волновым функциям от произведения определенного оператора на экспоненту действия [2, 3].
Для численных исследований векторных моделей элементарных частиц на решетке с существенным сокращением времени вычислений были предложены два алгоритма вычислений с генерацией случайных комплексных матричных полей, распределенных с заданным весом (плотностью вероятности): метод гибридного Монте-Карло [3, 4] и двухшаговый мультибозонный алгоритм [5, 6, 7]. Однако эти два альтернативных алгоритма, отличающиеся многими сложными процедурами и формулами, не были усовершенствованы оптимизацией параметров алгоритмов, не была доказана пригодность алгоритмов и не были сравнены их характеристики применительно к исследованию векторной U(1) модели элементарных частиц в разных ее фазах на решетке. Рассмотрение и решение этих проблем представляется весьма актуальным.
Целью данной работы являются анализ двух указанных алгоритмов применительно к векторной U(1) модели частиц на четырехмерной решетке пространства-времени, установление соотношений для оптимальных параметров и характеристик этих алгоритмов, вычисление с помощью каждого алгоритма значений трех важных корреляционных функций в Кулоновской фазе и фазе конфайнмента данной модели, выяснение пригодности рассматриваемых алгоритмов для данной модели по результатам вычислений, а также сравнение времени вычислений каждым алгоритмом.
Действие и корреляционные функции векторной U(1) модели
Согласно Вильсону [9], действие S[U, у, у ] векторной U(1) модели на четырехмерной решетке является суммой действия Sg[U калибровочного поля частиц - переносчиков взаимодействия и действия SF[U, у, у ] фермионных частиц (фермионов)
S[U, у, у ] = SG[U] + Sf[U, у, у ].
Здесь
Sg[U = PZ Re(1 - U„,), Sf[U, у, у ] =
x,|x,v
ц<у
Nf ____
TL'SM [ui, уy, (i)
f=1 x, y
U - решеточный компактный тензор на-
x,^v Г г
пряженностей калибровочного поля: U = U U+ U*. U ;
Ux л = exp(iAx л) - решеточное компактное калибровочное поле с вещественным потенциалом Axл в интервале (-п, п];
в = 1/e2 - обратный квадрат заряда; у^ и уХ - поля фермионов;
M[U] - фермионная матрица
M[U] = 5 - kV {(1 - у) U 5 л +
+ (1 + Y) U 5 л ;
4 1 л у,л y+л,*5
k = 1/(8 + 2m) - хоппинг-параметр;
m - фермионная масса;
Y - эрмитовы матрицы Дирака размером 4x4;
f - индекс поколения фермионов;
138
ЛЕСНОЙ ВЕСТНИК 6/2007
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Nf - число поколений фермионов; х, у - узлы четырехмерной конечной решетки пространства-времени с координатами nц = 0, 1, ..., N^ - 1; ц, v = 1, 2, 3, 4 - направления решетки; единичный вектор в положительном направлении;
N( - число узлов решетки вдоль направления ц.
Шаг решетки выбран равным a = 1.
В модели на решетке вычисление какой-либо корреляционной функции (O> выполняют методом континуального интеграла, состоящим в интегрировании соответствующей функции O[U] по всем полям U, у и у с весом exp(- S[U, у, у ]). При этом после интегрирования по полям у и у получается следующая формула [2-3, 9]
<O> = 1/Z O[U] exp(- SG [U])detN/M[U][dU], (2) где Z - нормировочная постоянная.
При исследовании моделей на решетке обычно вычисляют следующие важные корреляционные функции (O>: среднюю «энергию» калибровочного поля (Eg>, скалярный конденсат (у у> и «пионную» норму <П>. Эти величины рассчитывают по формуле (2), в которой зависимость O[U] полагают одной из следующих [1, 10]
Eg = 1/6V 2 Re0 - U,J. (3)
x,^,v
^<v
у у = 1/4V TrM-'[U], (4)
П = 1/4V Tr(y5M-1[U]y5 M-'[Uj), (5)
где V = N N2 N3 N4 - полное число узлов четырехмерной решетки пространства-времени, называемое объемом решетки; у5 = Y1 Y2 Y3 Y4 - киральная матрица Дирака; Tr - след в пространстве, являющемся прямым произведением пространства матриц Дирака и пространства узлов решетки x.
Метод гибридного Монте-Карло
Метод гибридного Монте-Карло HMC (hybrid Monte Carlo) предназначен для приближенных, но относительно быстрых вычислений корреляционных функций моделей на решетке с учетом детерминанта фермионной матрицы при четном числе поколений фермионных частиц N. В данном методе вводят вспомогательные поля P и х таким обра-
зом, что формула (2) переходит в следующее соотношение [4, 5]:
(O> = 1/Z O[U]exp(-H[U, P, x])[dUdP dx], (6) где H[U, P, x] - гамильтониан алгоритма, определенный формулой
H [U, P,x]= Sg [U]+12 p2 +
+f!x} (M [u~мри' x/
2
х,Ц
(7)
где Px ц, - вещественные компоненты поля P, называемого сопряженным импульсом калибровочного поля; f - поля с 4V комплексными составляющими, называемые псевдофермионами.
Поля U, P и x генерируют случайным образом с весом exp(-H[U, P, х]). Тогда формула (6) переходит в следующее соотношение для численных расчетов
O = lim — 2 O Ua) , (8)
\ f N^ Na=1 L J’ W
где O[U(a)] - значения функции O[U] при конфигурации поля U(a);
N - число конфигураций Ua\
Каждую конфигурацию полей U, P и х с указанным весом распределения получают по следующей схеме [4, 5]. Сначала генерируют случайную конфигурацию сопряженного импульса калибровочного поля P с весом
exp
2 P;,/ 2
х,Ц
и выбирают конфигурацию псевдофермионов X по формуле
X/=MU,
где П/ - случайный вектор с 4V комплексными составляющими, генерируемый с весом exp(-n/tnf); f = 1, ..., N/2.
Поля P и П/ генерируют методом стохастического вектора [3].
Затем выполняют процедуру молекулярной динамики, состоящую в получении промежуточной U' и вспомогательной P' конфигураций полей по известным конфигурациям U, P и х по следующим формулам Uk) = Uk-1) exp(/ Ат P(k-1/2) ),
х,ц х,ц А 4 х,Д/7
p(k+1/2) = p(k—1/2) + Ат F [ЦЮ]
где Ат - шаг «времени»;
N - число шагов «времени»; к = 1, ..., N; U0) = U ; UNt) = U' ;
7 7 т7 х,ц х,д7 х,ц х,д7
F [U] = -dH/dA .
х,ц L J х,ц
ЛЕСНОЙ ВЕСТНИК 6/2007
139
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Рассматриваемая схема заканчивается шагом Метрополиса, в котором полученную конфигурацию U принимают за последующую U с вероятностью wacc[U, U] = min{1, exp(H[U, P, X - H[P,P\X]}.
Если U не принято, то новая конфигурация совпадает с прежней U.
В процедуре молекулярной динамики и в шаге Метрополиса необходимо умножать псевдофермионы f на обратные матрицы M4 и [МД-1. Это умножение выполняют путем решения систем уравнений M2, = п с заданными п и неизвестными 2, векторами с 4V комплексными составляющими. Для ускорения решения таких систем используют «шахматное» разложение решетки на одинаковое число четных e (even) и нечетных o (odd) узлов с соответствующим разложением величин M, 2, и п
M =
1
M
MeC
1
п=
Ч ^ Uj
$=
Ч
Ч j
(9)
(10)
(11)
Тогда вместо указанной системы решают следующую эквивалентную ей систему для неизвестных векторов 2e и 2o с 2V составляющими:
Q = П -Moen ,
I ^e = П -Meo ^o ,
где Q - матрица размером 2V х 2V Q = 1 - M M .
Нахождение вектора 2o в системе (10) выполняют методом сопряженных градиентов [3]. Для этого верхнюю часть системы (10) путем умножения на матрицу Q' приводят к системе вида Ax = у, где A = Q'Q, x = 2o - неизвестный вектор, у - заданный вектор. Метод сопряженных градиентов CG (conjugate gradient) представляет собой следующую схему решения указанной системы путем построения векторов x быстро сходящихся к х:
X, = х + ag ; а = ||r ||2/g' Ag ; r = у - Ax ;
n+1 n n 11 n11 ° n °n7 n S n
gn+1 = rn+1 + Pngn; Pn = \\rn+l\\2/\\rЖ, (12)
где n = 0, 1, ...; ||a|| = Va'a; g0 = r0.
Выполнение этой схемы прекращают при достижении условия
4J < 5, (13)
где 5 - заданная величина.
5 = 5md в процедуре молекулярной динамики, 5 = 5acc в шаге Метрополиса, причем 5 << 5 << 1.
acc md
W < 2Pn/(1 +P2n)||r0||, р = (VC - 1)/ (VC + 1), С = Xmax/Xm,n. где X и X - максимальное и минималь-
max mm
ное собственные значения матрицы
A = QQ;
Z - ее число обусловленности.
Из указанного неравенства для нормы ||rn|| получим при Z >> 1 следующее соотношение для числа итераций NCG) схемы (12) метода CG до достижения условия (13):
NCG)in2. (14)
25
Двухшаговый мультибозонный алгоритм
Двухшаговый мультибозонный алгоритм TSMB (two-step multiboson algorithm) предназначен для приближенных быстрых вычислений в моделях на решетке при любом числе их поколений N. В данный алгоритм вводят вспомогательные поля мультибозоны Ф, мультибозонное действие SB^, U] и аппроксимирующие многочлены Pk(x) таким образом, что формула (2) преобразуется к виду [8]
(op ]det-1 P4 (Q [U ]Qp ]signdetf Qp ] (det-1 P4 @pl@ppgndetf Qp]
. (15)
Здесь обозначено:
<A>12 = 1/ZjA[U] exp(-SJU] -- Sb^, U]) det-1P2(Q'[U]Q[U])[dUdФ], (16) где SG[U] - действие по (1);
Q - матрица по (11);
sb№, u = £ф/ m [u] m [щф
j=1
' 1e Meo \
YsMoe (5 -Pj)le_ ;
Фj - поля мультибозоны с 4 V комплексными составляющими;
1 e, Meo и Moe- матрицы по (9); р. - комплексные корни многочлена
j P,(x2);
Pj(x), P2(x), P4(x) - положительные многочлены степеней n1, n2, n4 для грубой, промежуточной и достаточно точной аппроксимаций функции xNf2 на интервале x е [s, X]
Pj(x) « x-Nf/2, Pj(x)P2(x) = x-Nf/2,
Pl(x)P2(x)P4(x) = x-Nf2. (17)
140
ЛЕСНОЙ ВЕСТНИК 6/2007
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Интервал [в, к] содержит средние собственные значения (к . ) и (к ) матрицы Q}Q, усредненные по всем полям U с весом exp(-SG[UDdetwM[U].
Поля U и Ф генерируют случайным образом с весом exp(-SG[U] - £В[Ф, U]) х xdet-1/’2(Qt[U]Q[U]). Тогда формула (16) переходит в следующее соотношение для численных расчетов, аналогичное (8)
a = Nm -n 5 ACU“>]. (18)
Каждую конфигурацию полей U и Ф с указанным весом распределения получают путем выполнения двух шагов [6, 7, 8]. Для этого используют «шахматное» разложение решетки и действия £ДФ, U]) и SG[U] + ^В[Ф, U]) разбивают на слагаемые, содержащие составляющие полей Ф и U , и слагаемые без
j,x,r х,ц5
этих составляющих, гдеj = 1, ..., n1; x, д - узел и направление решетки; r - индекс строк или столбцов матриц Дирака.
В первом шаге каждую новую конфигурацию поля Ф' и промежуточную конфигурацию U получают из предыдущих конфигураций Ф и U путем выполнения следующих четырех процедур методами «тепловой бани» HB (heatbath) и верхней релаксации OR (overrelaxation) [6, 8]:
HB \
Ф' ,x,r =-Vj ,x,r/Aj ,x,r + S j,x,r,
OR \
Ф ' =-ф - 2V /A
j ,x,r j ,x,r j ,x,r / j ,x,r’
<r=KrFJ Fx,J)-2-
где V. = НТФ, U] не зависит от Ф ;
j,x,r L 7 J j,x,r
Fx^ = А[Ф, U] не зависит от U, д;
A не зависит от Ф и U;
j,x,r ’
E,jxr - комплексное случайное число, генерируемое методом стохастического вектора с весом exp (-A |S |2);
Пхд - вещественное случайное число в интервале (-п, п], генерируемое с весом exp(|F |cosn ).
Далее выполняют второй шаг алгоритма TSMB - метод Метрополиса [7, 8]. В этом шаге полученную конфигурацию поля ^/принимают за последующую U с вероятностью
waJU, U] = min(1, detP2 (Q^[U']Q[U'])/ /det-1P2 (Q\U]Q[U\). (19)
Если U не принята, то новая конфигурация совпадает с прежней U. Здесь отношение детерминантов вычисляют по формуле [8] det-1P2 (Q\U']Q[U'])/detlP2 (Q^[U]Q[U]) =
KexpK^QWMU])^ + Чп))л, (20) где
(B)n = 1/ZjB[n] exp(-nfn[dn]; (21)
U - вектор с 2V комплексными составляющими;
S = PQUQmn;
P3(x) - положительный многочлен степени n аппроксимирующий с высокой точностью функцию P2-1/2(x) на интервале x е [в, к]
P3(x) = P2-1/2(x). (22)
При вычислениях величин (O) по (15) с (16) и (18) используют следующую формулу [8]
det-1p4(QtQ) = <exP(-nP4(QtQ) П + (23)
где усреднение выполнено согласно (21).
Коэффициенты аппроксимирующих многочленов Pk(x) с заданными степенями nk, где к = 1, 2, 4, на интервале x е [в, к] вычисляют интегральным методом наименьших квадратов. Коэффициенты многочлена P3(x) с заданной степенью n3 находят итерационным методом, аналогичным методу Ньютона касательных. Вычисление матриц Pk(QQ), где к = 2, 3, 4, в формулах (20) и (23) осуществляют по следующим соотношениям, в которых x заменяют на Q}Q
Пк
Pk(x)=Zcf °)(x) , Oj+1(k) =
j=0
= (x + Pj(k)) Oj(k) (x) + yj-1(k) °j-1(k)(x), где Oj(k) (x) - ортогональные многочлены степени j.
Собственные значения к и к
min max
матрицы A = QQ, необходимые для выбора границ интервала [в, к], для каждой конфигурации поля U находят явно повторяемым методом Ланцоша. Этот метод представляет собой итерации, включающие:
1) для итерации с номером n последовательное вычисление M единичных взаимно ортогональных векторов e(n\ где j = 0, 1, ..., M - 1, по формуле
Ae(n) = В(n) e.,,
j w j+1
(n) + a(n) e(n) + В (n) e (n)' j j *7-1
j-1
2) вычисление либо к (n), либо к
min m
(n)
собственного значения трехдиагональной
ЛЕСНОЙ ВЕСТНИК 6/2007
141
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
матрицы Т(п) размером M х M с элементами а(п) на главной и в-п) на соседних диагоналях, а также вычисление соответствующего единичного собственного вектора Уп) с составляющими Ап);
3) нахождение для следующей итерации вектора
ет)= i1...
J=0 j j
Средняя вероятность (wacc) принятия новых полей U равна отношению числа конфигураций, принятых за новые, к полному числу генерируемых конфигураций. Учитывая соотношения (19) и (20) получим для алгоритма TSMB при ((AE)2) < 1:
<w->=1 •
где
<(AE)2) ~ F|P2(x) - 1||2,
(24)
(25)
AE = ^ши-миге - пЧ S определено в (20);
(C) = 1/Z\Cpi2exp(-SG[U] - 5ДФ, U] -- пЧ [dUdФ dUWdn];
С = С[иФ'иФц]; p12 = p12[UФ'иФ] -плотность вероятности перехода от конфигурации полей U и Ф к последующей конфигурации U и Ф .
Параметры метода гибридного Монте-Карло
Точность и продолжительность расчетов методом HMC зависят от значений его следующих исходных параметров: Ат, A, 5md
и 5
acc
Для получения необходимых соотношений рассмотрим следующую матрицу Ц размером 4V х 4V
Ц=
fd2 H Y fd2 H Л
nV2\
ydA2 j
ydA2 j
где H - гамильтониан по (7);
A - потенциал поля U; усреднение выполнено, как в (6).
Оценим норму этой матрицы ||Ц|| при Nf = 2. Для этого используем свойства симметрии U(1) модели относительно дискретных сдвигов и поворотов, формулы
(ХХТ)Л = MMи XjMM) = X2[(MM)1/2], а также сумму по j заменим интегрированием. В результате получим следующее соотношение
|Ц| ~ hr M tM )-
Ч ~ V f А. (26)
A kj/ АХX
где X, X и X - собственные значения мату mm max
рицы (MM)1/2;
Z - число обусловленности Q}Q
Z = <Xmax(QtQ))/<XmIn(QtQ)). (27)
При исследовании модели свободного бозонного поля методом HMC ранее были получены следующие соотношения для обеспечения достаточно точных и быстрых вычислений
V| Ц |2 (Ат)4 <1, NtAt ~1, (28)
где
Ц
д2 н
~д£
- матрица, не зависящая от бозонного поля A;
H - гамильтониан бозонного поля; V -объем решетки.
Мы применим к U(1) модели соотношения (28), в которых положим Ц= Ц и учтем (26). Тогда получим следующие соотношения для параметров Ат и А
Ат<-
1
NT >V1/4.
(29)
V с)'4’
Для остальных параметров метода HMC принимаем следующие соотношения
5 ~ 1/V, 5 ~ 1/V2. (30)
Из (14) и (30) получим следующее соотношение для среднего числа итераций (ACG)) в методе CG применительно к методу HMC
<ACG)) ~ VZlnV. (31)
Параметры двухшагового мультибозонного алгоритма
Для расчетов алгоритмом TSMB должны быть заданы его следующие исходные параметры: границы s и X интервала аппроксимаций функции x~Njf!2 по (17) и степени п п п3 и п4 аппроксимирующих многочленов.
Границы s и X выберем следующим образом
s = 0,5<Xmm); X = (1,2 - М) <Xmаx), (32)
где (Xmin) и (Xmax) - средние собственные значения матрицы QQ, усредненные, как в (15).
142
ЛЕСНОЙ ВЕСТНИК 6/2007
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Такой выбор учитывает необходимость нахождения значений (A . > - oA и
(kmax> + 0Amax внутри интервала k AL где 0Amin
и oAmax - среднеквадратичные отклонения.
Достаточно точные и быстрые вычисления алгоритмом TSMB обеспечиваются при условиях
<0 ~ 1, (det-P (QQ)> ~ 1.
При этих условиях из (23) и (24), а также из условий достаточно точных аппроксимаций (17) и (22) получим следующие соотношения
_1_ W ’
X
N
2P1 (X)P2 (X)-1
-, (33)
V
где |]fx)|| определено по (25).
Для оценки значений степеней nk многочленов Pk(x) введем вспомогательные многочлены Pk (х) степеней Пк, составленные из многочленов Чебышева Тп +1 (х), со следующими аппроксимациями:
P1 (X)~ P1 (X)~ xNf!2 , P2 (X)~ P1 (X)P2 (X)~ xNf!2 ,
P (х> P2 (x)P (х)=1,
P4 (X)~ P1 (X)P2 (X)P4 (X)= хЫ^ , (34)
где
Pk (X)= ^ X
T
1-
2 x-(A'+s') A'-s'
T +1-
А'+£Л
nk+11 A'-s',
Tk+1 (X)= cos[(nk +1>rccosX] . (35)
Согласно нашим расчетам при Nf = 2, многочлены P1 (х) и P1(x) практически одинаково аппроксимируют функцию x~Nf2 при выполнении следующих условий
nl=n,, s'-VsA, А' = А . (36)
Из (27) и (32) следует соотношение в
виде
Nf/2
A/s ~ Z. (37)
Из (33, 34, 35-37) получим
-2(ч+1)-1/4 ~NL
P (х)-1
( A'+s' ^ I A'-s'
откуда имеем соотношение для степени n1
n1 ~ Z1/4lnV. (38)
Наши расчеты при Nf = 2 показали, что многочлены P2 (х) и P1(x)P2(x) практически одинаково аппроксимируют функцию xNf2 при условиях
П2 = ц+ n2, s'-s, A'=A. (39)
Из (33, 34, 35, 37, 39) получим
х
Nf12
P (х)-1
' I
A'+s
n2+1l A'-
-2(1+n, +1) 12 1
V
откуда имеем соотношение для степени п2
e
n2 ~ Z1/2lnV. (40)
Результаты аналогичных расчетов значений многочленов P3 (х), P4 (X), P2(X)P32(X) и P1(x)P2(x)P4(x) с аппроксимациями (34) и последующий анализ приводят к следующим соотношениям для степеней п3 и п4
П > ^ П4 > П2. (41)
Производительности алгоритмов
Производительностью алгоритма вычислений конкретной корреляционной функции (O> модели на решетке называем величину P, равную
P=(42) где (Nop> - среднее число компьютерных операций, необходимых для перехода от одной конфигурации полей к последующей;
(т .nt> - интегральное время автокорреляций, равное среднему расстоянию между ближайшими статистически независимыми значениями 0[^(а)]этой корреляционной функции [3].
Это время вычисляют несколькими способами, например, способом «суммирования» (ты) рассчитывают по формуле
где
т \=1+1 £ МЙ
int/ 2 W £1 Г0 (0)
Г (а)=—1— х oW N-а
(43)
N-а( 1 N-а V 1 N
х ЦО,- — ZOk ON— £ Ok |;
j=1 \ N-а k=1 Л N-аk=а+l
1<< W<<N; Оа = O[№4; остальные обозначения, как в (8).
Время вычисления t величины (O> с заданной точностью, практически равное времени генерации N конфигураций полей Ц(а), определяется формулой
t = <WN = WP (44)
где (tCPU> = tc(Nop> - среднее компьютерное время перехода от одной конфигурации полей к последующей;
ЛЕСНОЙ ВЕСТНИК 6/2007
143
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
tc - среднее время выполнения одной компьютерной операции;
N = N/(t t) - число статистически независимых значений O .
а
Лучшим алгоритмом, т.е. алгоритмом с меньшим временем вычислений (О) при одинаковых точности вычислений, параметрах модели N, V, в, k и числе N , является алгоритм с большей производительностью P.
Получим соотношения для производительностей рассматриваемых алгоритмов. В случае метода HMC на модели свободного бозонного поля на решетке было получено следующее соотношение
(xmtHMC) ~ (NAt)-2. (45)
Мы применим к 7(1) модели это соотношение, в которое подставим параметры по (29).
В случае алгоритма TSMB справедливо следующее соотношение
<т JSMB ) ~ (min|Imp 2|)-1, (46)
где р. определены в матрице M к (16). Заменим, согласно (34), многочлен P1(x) на р (x) и с учетом (35) получим соотношение для корней р. многочлена р (x). Подставляя это соотношение в (46) и учитывая (27, 32, 36, 37), получим
(
(ТГ) ~
/s'sin
2п
V1
А (А
Us'
aV2
А +1
( л V/2 '
А X
14
(47)
Для каждого алгоритма вычислений в 7(1) модели имеем
(NopHMC) = 4V7(1 + 2(N-CG))) +
ОР + 8V(Ncg)) ~ V(Ncg))Nt, (48)
(NJsmb) = 2V(128a1 + n2 + n)~ Vn2. (49) Из (42) с учетом (29), (31), (45) и (48) для метода HMC, и с учетом (38, 40, 47, 49) для алгоритма TSMB получим следующие соотношения для производительностей этих алгоритмов применительно к 7(1) модели:
Phmc ~ ртшв ~ U^W. (50)
Отсюда для отношения производительностей имеем
P /P ~ V1/4/lnV (51)
Результаты вычислений и сравнение алгоритмов
Нами были разработаны, отлажены и применены для расчетов программы вычисле-
ний средних значений корреляционных функций рассматриваемой векторной 7(1) модели на четырехмерной решетке как методом HMC, так и алгоритмом TSMB. По этим программам были вычислены указанные в п. 2 величины (Eg), (7v) и (П) данной 7(1) модели в областях Кулоновской фазы и фазы конфайнмента. Вычисления проведены на суперкомпьютере VPP-300 (Ганновер, Германия).
Были заданы следующие параметры 7(1) модели, одинаковые при вычислениях каждым алгоритмом: число поколений фермионов Nf = 2; объем решетки V = 63 х 12 при N = N2 = N3 = 6 и N3 = 12; параметры в = 2 и k = 0,130 для Кулоновской фазы, в = 0 и k = 0,238 для фазы конфайнмента; число генерируемых полей Uf-a)N = 10000; периодические граничные условия по направлениям д =1, 2, 3 и антипериодические по направлению д = 4 для фермионных полей, а также периодические условия по всем направлениям д = 1, 2, 3, 4 для калибровочных полей. Выбранные значения в и k находятся вблизи линии раздела фаз, что необходимо для правильного описания свойств частиц. Кроме того, интересно узнать характеристики алгоритмов при «трудном» для них значении в = 0 [10].
Для расчета значений исходных параметров каждого алгоритма в каждой фазе сначала были вычислены предварительные значения (Xmin) и (Xmax) матрицы Q}Q путем усреднения X . и X по полям 7 с весом exp(-SG[7]). При этом использован явно повторяемый метод Ланцоша, а генерация конфигураций 7 выполнена методом HB при Ф = 0. По формулам (27, 29, 30) с этими значениями <Xmin) и (Xmax) были найдены предварительные значения параметров метода HMC. Затем методом HMC с этими предварительными параметрами была снова выполнена генерация полей 7, но с весом exp^S^T^detTMCT. C помощью этих полей были получены следующие окончательные значения, одинаковые в каждом алгоритме: (Xmin) = 0,13(1) и (Xmax) = 1,63(1) в Кулоновской фазе, и 0,0005(1) и 6,59(1) в фазе конфайнмента.
Используя эти окончательные значения <Xmin) и <Xmax), по (27, 29, 30) для метода HMC, и по (27, 32, 38, 40, 41) для алгоритма TSMB были рассчитаны значения парамет-
144
ЛЕСНОЙ ВЕСТНИК 6/2007
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ров алгоритмов (табл. 1). В случае алгоритма TSMB интегральным методом наименьших квадратов найдены коэффициенты аппроксимирующих многочленов Pk(x), к = 1, 2, 3, 4, а также найдены корни многочлена Pj(x2).
Каждым из расматриваемых алгоритмов с использованием формул (3-5, 8, 15, 16, 18, 23) были рассчитаны корреляционные функции (Eg), ( у V) и (П) данной U(1) модели в двух ее фазах. Кроме того, были рассчитаны значения величин (wacc), (хт) по (43), (tCpu), а также отношение производительностей по формуле, следующей из (42) с учетом (44)
P P = (t HMC)(T HMC)/(t TSMB)(T TSMB) TSMB HMC \ CPU int /'\lCPU int /■
Результаты выполненных расчетов приведены в табл. 2. Значения каждой из трех корреляционных функций в каждой из двух фаз модели, полученные с применением двух разных алгоритмов, оказались одинаковыми. Одинаковой оказалась также точность вычислений. Это свидетельствует о пригодности каждого из рассматриваемых алгоритмов, отличающихся многими сложными процедурами и формулами, для вычислений с существенным сокращением времени в U(1) модели по Вильсону на четырехмерной решетке. Кроме того, это подтверждает правильность установленных соотношений для исходных параметров каждого алгоритма. Особо отметим совпадение результатов вычислений каждым алгоритмом и пригодность этих алгоритмов для «трудного» случая в = 0 в фазе конфайнмента.
Полученные средние значения (wacc) для каждого алгоритма находятся в пределах 0,48-0,94, что соответствует принятой оценке (wacc) ~ 1. Значения (тп) для каждой корреляционной функции оказались одного порядка величины для разных алгоритмов, причем для величин (Eg) и ( у у) в фазе конфайнмента при в = 0 они оказались на полтора порядка больше, чем в Кулоновской фазе при в = 2, а для корреляционной функции (П) - одного порядка величины. Значения (tCPU) оказались одного порядка величины в каждой фазе при разных алгоритмах, но в фазе конфайнмента они были в 5-8 раз больше, чем в Кулоновской фазе.
Отношения производительностей алгоритмов для шести рассмотренных случаев (три корреляционные функции и две фазы)
оказались в пределах PTSMBPHMC = 0,5-1,8. Этот диапазон соответствует аналитическому соотношению (51), по которой для данного объема решетки V = 63 х 12 отношение PTSMBl PHMC ~ 0,91. При увеличении объема решетки производительность каждого алгоритма, согласно (50), уменьшается, а время вычислений, соответственно, увеличивается. Но отношение PTSMB/PHMC, согласно (51), возрастает, и алгоритм TSMB становится лучшим по производительности, т.е. с меньшим временем вычислений, чем метод HMC.
Т а б л и ц а 1
Значения исходных параметров двух алгоритмов для двух фаз U(1) модели
Фаза HMC
Ат &md &acc
Кулоновская 0,025 40 10-3 10-7
конфайнмента 0,01 10 10-3 10-7
фаза TSMB
S X «1 «2 «3 «4
Кулоновская 0,025 2,5 6 36 48 200
конфайнмента 0,000225 9 50 360 450 500
Т а б л и ц а 2
Значения корреляционных функций
(O) в двух фазах U(1) модели, вычисленные двумя алгоритмами, и характеристики алгоритмов
Параметры (Eg) (vv) (П)
Кулоновская фаза
(OHMC) 0,1332(1) 0,9381(1) 1,378(1)
(QTSMB) 0,1331(1) 0,9379(1) 1,376(1)
(w HMC) acc 0,94(1)
(WaccTSMB) 0,48(1)
КГ) 3,2(3) 2,0(2) 25(4)
(Ti„,TSMB) 3,0(3) 2,8(2) 50(8)
(tCPU^ с 15,1(2)
/t tsmb\ n VCPU /> ^ 8,96(2)
P TSMBlPHMC 1,8(2) 1,2(2) 0,8(2)
фаза конфайнмента
(OHMC) 0,939(1) 0,95(1) 13,9(2)
(QTSMB) 0,938(1) 0,96(1) 13,7(2)
(w HMC) acc 0,72(1)
(WaccTSMB) 0,68(1)
KD 65(7) 60(7) 35(5)
(T JSMB) 120(20) 125(15) 45(5)
p 1 о 76(1)
Zt tsmb\ n VCPU Л ^ 69(1)
P TSMB1PHMC 0,6(1) 0,5(1) 0,9(1)
ЛЕСНОЙ ВЕСТНИК 6/2007
145
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Рассмотрены и реализованы два алгоритма для векторной 6(1) модели на четырехмерной решетке пространства-времени: метод гибридного Монте-Карло и двухшаговый мультибозонный алгоритм. В каждый алгоритм входят генерация калибровочных полей переносчиков взаимодействия при одинаковом весе распределения с детерминантом фермионной матрицы, метод стохастического вектора и «шахматное» разложение решетки.
Данные алгоритмы отличаются многими входящими в их состав формулами, методами, процедурами и подходами. В метод гибридного Монте-Карло входят использование вспомогательных полей псевдофермионов и сопряженного импульса калибровочного поля, процедура молекулярной динамики для обновления значений полей, метод сопряженных градиентов и шаг Метрополиса для гамильтониана алгоритма. Двухшаговый мультибозонный алгоритм включает использование вспомогательных полей мультибозонов, методы «тепловой бани» и верхней релаксации для обновления полей, применение многочленов для аппроксимации детерминанта фермионной матрицы, явно повторяемый метод Ланцоша для собственных значений этой матрицы и шаг Метрополиса.
Установлены соотношения для оптимальных параметров и характеристик этих алгоритмов применительно к рассматриваемой 6(1) модели. Это зависимости для шага «времени», числа шагов «времени» и параметров остановки метода сопряженых градиентов в методе гибридного Монте-Карло, соотношения для границ интервала аппроксимаций и для степеней четырех аппроксимирующих многочленов в случае двухшагового мультибозонного алгоритма, а также зависимости для производительности каждого алгоритма.
Разработаны, отлажены и реализованы программы вычислений корреляционных функций 6(1) модели каждым из указаных алгоритмов. С помощью каждого алгоритма в областях Кулоновской фазы при параметре модели в = 2 и фазы конфайнмента при «трудном» значении в = 0 вычислены значения следующих важных корреляционных функций: средней «энергии» калибровочного поля, скалярного конденсата и «пионной» нормы. Значения корреляционных
функций в каждой фазе, рассчитанные разными алгоритмами, оказались одинаковыми при одинаковой точности.
Установлено аналитически и подтверждено численными расчетами, что для выбранной решетки пространства-времени с числом ее узлов 63 х 12 время вычислений корреляционных функций с помощью каждого алгоритма является примерно одинаковым. При увеличении объема решетки время вычислений каждым алгоритмом возрастает, но алгоритмом с меньшим временем вычислений становится двухшаговый мультибозонный алгоритм.
Указанные алгоритмы с параметрами по установленным зависимостям целесообразно использовать при исследованиях моделей элементарных частиц методом решетки.
Автор выражает благодарность за плодотворные дискуссии д-ру физ.-мат. наук И.Л. Боголюбскому, проф. И. Монтвею, канд. физ.-мат. наук В.К. Митрюшкину и проф. М. Мюллеру-Пройскеру. Работа выполнена при финансовой поддержке грантом Правительства Москвы и Департамента образования города Москвы в области наук и технологий в сфере образования за 2003 и 2005 гг.
Библиографический список
1. Кройц, М. Кварки, глюоны и решетки / М. Кройц.
- М.: Мир, 1988.
2. Славнов, А.А. Введение в квантовую теорию калибровочных полей / А.А. Славнов, Л.Д. Фаддеев.
- М.: Наука, 1988. - 272 с.
3. Montvay I., Monster G. Quantum Fields on a Lattice. Cambridge: University Press, 1994.
4. Duane S., Kennedy A., Pendleton B., Roweth D. Hybrid Monte Carlo // Phys. Lett. B. 1987. V 195. P. 216-222.
5. Gottlieb S., Liu W., Toussaint D., Renken R., Sugar
R. Hybrid molecular dynamics algorithms for the numerical simulation of quantum chromodynamics // Phys. Rev. D. 1987. V 35. P. 2531-2542.
6. Loscher M. A new approach to the problem of dynamical quarks in numerical simulations of lattice QCD // Nucl. Phys. B. 1994. V 418. P. 637-648.
7. Borici A., de Forcrand Ph. Systematic errors of Loscher’s fermion method and its extensions // Nucl. Phys. B. 1995. V. 454. P. 645-662.
8. Montvay I. An algorithm for gluinos on the lattice // Nucl. Phys. B. 1996. V 466. P. 259-284.
9. Wilson K.G. Confinement of quarks // Phys. Rev. D. 1974. V. 10. P. 2445-2459.
10. Hoferichter A., Mitrjushkin V.K., Mьller-Preussker M., Stoben H. Dynamical Wilson fermions and the problem of the chiral limit in compact lattice QED // Phys. Rev. D. 1998. V 58. P. 114505-114510.
146
ЛЕСНОЙ ВЕСТНИК 6/2007
