CC BY
28
Математическое моделирование. Оптимальное управление Вестник Нижегородского университета им. Н.И. Ло(Зачевского, 2011, № 3 (2), с. 62-66
УДК 519.6+530.182+534.1
АНАЛИЗ ТОЧНЫХ ПЕРИОДИЧЕСКИХ РЕШЕНИЙ В БОЛЬШИХ ОСЦИЛЛЯТОРНЫХ ЦЕПОЧКАХ С ИСПОЛЬЗОВАНИЕМ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ
© 2011 г. М.В. Иванченко
Нижегородский госуниверситет им. Н.И. Лобачевского
ivanchenko@rf.unn.ru
Поступила в редакцию 28.10.2010
Разработан алгоритм численного нахождения и анализа линейной устойчивости точных периодических решений в цепочках нелинейных осцилляторов на базе параллельных вычислений. Достигнутое ускорение позволило впервые получить такие решения в цепочках длиной более тысячи элементов, показать их локализацию в пространстве нормальных мод системы, исследовать устойчивость, сопоставить результаты с теоретическими предсказаниями. Обсуждается прикладное значение полученных результатов.
Ключевые слова: нелинейная динамика, параллельные вычисления, динамика ансамблей.
Введение
Исследование динамических возбужденных состояний является одной из фундаментальных задач физики твердого тела. Основная трудность в классическом многочастичном описании нелинейных систем заключается в установлении связи между возбужденными состояниями и характерными траекториями динамической системы в фазовом пространстве. Поскольку в общем случае гамильтоновы системы неинтегрируемы, динамика будет хаотической практически для любых начальных условий. Это, в свою очередь, означает релаксацию возбужденных состояний на достаточно больших временных масштабах, поскольку последние описываются регулярными колебаниями. Однако Колмогоров, Арнольд и Мозер показали, что регулярные колебания на инвариантных Ы-мерных торах могут сохраняться и в неинтегри-руемых системах с конечным числом степеней свободы [1]. Особенно важную роль играют простейшие из низкоразмерных инвариантных структур - периодические траектории. Они существуют даже в сильно хаотических системах и плотно заполняют хаотическую часть фазового пространства несмотря на их нулевую меру.
Периодические траектории являются базовым классом регулярных траекторий, описывающих возбужденные состояния. Одним из примеров таких состояний являются дискретные бризеры -периодические траектории в нелинейных колебательных решетках, локализованные в прямом
пространстве [2]. Их аналог - периодические траектории, локализованные в пространстве нормальных мод, - д-бризеры, которые обеспечивают баллистический транспорт энергии и аномальную теплопроводность в низкоразмерных решеточных структурах [3]. Известны численные методы, позволяющие найти эти периодические решения с очень высокой точностью [2,4]. Они ярко иллюстрируют роль численных исследований в современных вопросах физики нелинейных возбужденных состояний.
Основной проблемой численного анализа является возрастание трудоемкости вычислений с увеличением размеров решетки. Так, до настоящего времени д-бризерные решения были получены для цепочек длиной не более нескольких сотен элементов. Анализ их устойчивости является еще более трудоемкой задачей. Размер реальных физических систем гораздо больше; например, нанотрубки имеют порядка десятков тысяч атомов в длину. С необходимостью исследования моделей, близких к реальным, и связана актуальность разработки численных методов отыскания периодических траекторий и исследования их устойчивости с использованием параллельного программирования.
Математическая модель
Разработанный численный алгоритм решает задачу продолжения точных периодических траекторий линейной системы в нелинейную область для произвольных потенциалов взаимо-
действия в цепочке и собственных колебаний осцилляторов, а также не содержит ограничений на вид (свойства) этих периодических траекторий. В данной работе мы рассматриваем простейшую динамическую модель атомарной цепочки типа Ферми-Паста-Улама с произвольным порядком нелинейности в потенциале взаимодействия и уравнениях движения [5]:
х = х , — 2х + х , + х(х , — х У +
п п+1 п п-1 ' П+1 я/ '
+ Х( хя-1 — хп)
У-1
(1)
где хп - координата п-го осциллятора, х - коэффициент нелинейности, у - показатель нелинейности, граничные условия выбраны фиксированными: х0 = хм+1 = 0. Каноническое преоб-2 N
разование хп (?) = '
71 в, У)
пап
8ш--------- дает
N +1 N +1
решение линейной задачи (% = 0) в амплитудах нормальных мод ^, частоты которых равны
В нелинейной системе такое
Ю„ = 281И------
4 2(Ж +1)
преобразование приводит к динамическим уравнениям для ^-осцилляторов с нелинейной связью:
X®,
Qq + ®д -
N
х 2 ....?т-1 П ®9к О9к ’
Ї1 .-.?у-1=1 к=і
[2( N + 1)]у/2-1
у-1
менить теорему Ляпунова [6], которая гарантирует возможность продолжения периодической траектории в нелинейную область при условии фиксированной энергии.
Реализация численного метода продолжения д-бризеров основывается на том, что периодическая траектория отвечает неподвижной точке в отображении, задаваемом секущей Пуанкаре в многомерном фазовом пространстве. Отображение может быть сконструировано следующим образом. В качестве секущей Пуанкаре
выбирается Х[2(N+1)/?0] = О, Х[N+1)/го ] > 0 ; начальные условия, отвечающие возбуждению в виде линейной моды q0 задаются на этой секущей. Система дифференциальных уравнений численно интегрируется до повторного пересечения фазовой траекторией секущей Пуанкаре, которое задает отображение начальных импульсов и координат на себя: ум = р(уп). Неподвижная точка отображения является решением системы нелинейных уравнений О(у) = = Р(у) - у = 0, поиск которого осуществляется обобщенным методом Ньютона-Рафсона. Для этого параллельно с интегрированием исходной системы (1) производится интегрирование линеаризованных уравнений (1) для 2Ы линейно независимых единичных векторов возмущений, и разность их конечных и начальных значений
дб (у)
формирует матрицу Ньютона: N(у) = -
(2)
ду
где коэффициенты Бя определяют взаимодействие между конкретными модами [5]. Нормальные моды не являются точными периодическими решениями нелинейной задачи, однако можно ожидать, что они могут быть продолжены в нелинейную область, по крайней мере, при малости коэффициента нелинейности.
Последовательный алгоритм
Последовательный алгоритм отыскания периодических траекторий в нелинейной цепочке состоит в продолжении нормальных мод линейной задачи в нелинейный режим [4]. Исходное возбуждение задается в виде моды Qq с
энергией Еч , что отвечает периодической траектории в многомерном фазовом пространстве линейной системы. Линейный спектр системы (1) удовлетворяет критерию нерезонансности ю Ф пю для всех а ф ап. Это позволяет при-
Я Яо II0
Решая систему линейных уравнений ^ (у) = = N(у)(у - у') методом Гаусса, находим новую оценку координат неподвижной точки у. Итерации повторяются до тех пор, пока не будет достигнута заданная точность ||р(у) - У /IУ «=•
Заметим, что если начальные условия находятся на периодической траектории, то матрица Ньютона N(у) будет являться матрицей Флоке для этой траектории, собственные числа которой представляют собой мультипликаторы и определяют ее линейную устойчивость. В нашем случае мультипликаторы вычисляются как собственные числа матрицы Ньютона, полученной на последней итерации обобщенного метода Ньютона-Рафсона.
Изложенный метод был успешно применен к анализу зависимости локализации и устойчивости д-бризеров от порядка нелинейности в уравнении (1) [5], однако вычислительная трудоемкость, возрастающая с увеличением размера системы, ограничивала доступный размер цепочки несколькими сотнями осцилляторов.
Параллельный алгоритм
Оценим вычислительную трудоемкость последовательного алгоритма для цепочки длиной N. На практике сходимость вычислительного метода начинает сильно зависеть от размера системы только при относительно больших значениях нелинейности/энергии системы, примерно соответствующих делокализации q-бризера в модовом пространстве. Тогда число итераций может существенно возрастать, либо сходимость исчезает. Мы будем работать в режиме хорошей локализации и быстрой сходимости, что на практике означает не более 5-10 итераций для достижения относительной ошибки в определении координат неподвижной точки отображения Пуанкаре менее 10-7.
Каждая итерация состоит из двух основных этапов: интегрирование системы (1) и 2N N-мерных векторов возмущений на периоде времени между последовательными пересечениями и решение линейной системы уравнений методом Гаусса. Первый этап имеет сложность O(N3), поскольку размерность интегрируемой полной системы (1) и линеаризованной системы возмущений равна 2N(2N+1), а время интегрирования растет как период линейной моды q0: Х. Второй этап также имеет сложность O(N3). В целях ускорения вычислений будем масштабировать номер зародышевой моды q0 с увеличением N: q0 / N = const, так что трудоемкость
первого этапа сокращается до O(N2), помня, вместе с тем, о фактической сложности исходной задачи.
Все численные эксперименты проводились на кластере ННГУ общей производительностью 2.7 терафлопс в составе 64 двухпроцессорных серверов на основе процессоров Intel Xeon 5150 2.66 GHz, 4GB DDR2-667 FB-DIMM RAM, 80 GB HDD, 2-х Gigabit Ethernet. Результаты сопоставления вычислительных затрат первого и второго этапов итерации последовательного алгоритма (Tmm, Tg) представлены на рис. 1. Они
согласуются с теоретическими оценками и показывают, что основные вычислительные затраты приходятся на интегрирование исходной и линеаризованных систем в исследовавшемся диапазоне N. Это обусловило реализацию параллельного алгоритма только для первого этапа итерации. Параллельная реализация метода Гаусса является решенной задачей [7] и может быть также использована при необходимости.
Первый этап итерации обобщенного метода Ньютона-Рафсона является очень удобной для распараллеливания задачей, поскольку состоит
Рис. 1. Зависимость времени выполнения этапов алгоритма от размерности задачи и числа процессов
в интегрировании исходной системы и 2Ы независимо эволюционирующих возмущений. Сбалансированная нагрузка и полное исключение затрат на коммуникацию между р процессами достигается разбиением задачи на р независимых частей, каждая из которых представляет собой интегрирование системы (1) в совокупности с порядка [2Ы/р] линеаризованных систем размерностью 2Ы каждая. После завершения процедуры интегрирования конечные состояния векторов возмущений собираются в основном процессе, где в рамках последовательной части алгоритма вычисляется новая точка на секущей Пуанкаре. Таким образом, вычислительная сложность численного метода с применением параллельного программирования будет порядка 4Ы2/р+аЫ, где а<<1, поэтому ускорение может быть весьма близко к р, а эффективность - к 1.
Результаты
Рассмотрим сначала зависимость времени, затрачиваемого на вычисление матрицы Ньютона и решение линейной системы, от размеров системы для параллельного алгоритма, выполняемого на 64 процессорах (рис. 1). Видно, что этап интегрирования существенно ускоряется и приближается по порядку малости к времени решения линейной системы. Последнее слегка увеличивается за счет затрат на сборку матрицы Ньютона и рассылку нового значения стартовой точки другим процессам для следующей итерации метода.
Зафиксируем теперь размер системы N = 256 и qo = 24 и будем исследовать зависимость ускорений для одного этапа итерации Бтт =
= Ттт (1)/ Ттт (р) и полного времени работы ал-
Рис. 2. Зависимость ускорения и эффективности параллельного алгоритма для времени вычисления матрицы Ньютона и полного времени решения задачи от числа процессов
Рис. 3. Распределение энергии q-бризеров в пространстве нормальных мод в зависимости от размеров системы и номера центральной моды. Сплошными линиями показана теоретическая оценка (3)
горитма £у = Гу (1)/ Гу (р), а также соответствующей эффективности Етт,у = Бтт,у / р [7].
Как следует из результатов, представленных на рис. 2, ускорение близко к максимальному, а эффективность - к 1 для сравнительно небольшого числа процессов. Для числа процессов порядка 100 потеря эффективности становится заметной. Сравнение временных затрат на этапы интегрирования и решения линейной системы показывает, что этот эффект вызван уменьшением первых вследствие распараллеливания.
Обсудим в заключение физическое содержание численных результатов. Впервые удалось найти точные периодические решения, локализованные в пространстве линейных мод - д-бризеры - в нелинейных цепочках размером более тысячи осцилляторов. Полученные профили распределения энергии в модовом пространстве для у = 4 (рис. 3) согласуются с теоретически предсказанной зависимостью [4]
Е,
(2 п+1) q^
= Е
qo
2п
А3ХЕ (N +1)Л
2 2 П %0
(3)
Получены численные результаты по устойчивости данных решений, которые свидетельствуют о линейной неустойчивости для N = = 2048, 4096, что согласуется с теоретическим результатом для порога неустойчивости 6%Едо (N +1) / п2 = 1 [4], что при данных значениях параметров дает N > 1644.
Заключение
В работе описан последовательный алгоритм, позволяющий находить точные периодические траектории в цепочках нелинейных осцилляторов общего вида путем продолжения известных решений для линейного случая в нелинейную область и проводить численный анализ их устойчивости. Разработан параллельный алгоритм для решения этой задачи. Исследованы ускорение и эффективность, даваемые этим алгоритмом. Использование параллельного алгоритма позволяет находить точные периодические решения в цепочках существенно больших размеров, чем было возможно ранее. Работа алгоритма продемонстрирована на примере задачи о д-бризерах, периодических траекториях, экспоненциально локализованных в пространстве мод и играющих принципиальную роль в физических процессах термализации и переноса энергии в кристаллических системах. Полученные профили распределения энергии и характеристики устойчивости согласуются с теоретическими результатами. Разработанный алгоритм открывает новые возможности численного исследования периодических траекторий нелинейных цепочечных систем и будет использован для решения прикладных физических задач.
Автор благодарит С. Флаха, О.И. Канакова и К.Г. Мишагина за плодотворные дискуссии.
Список литературы
1. Arnold V.I. Mathematical Methods of Classical Mechanics. New York: Springer-Verlag, 1989.
2. Flach S. and Willis C.R. Discrete breathers // Phys. Rep. 1998. V. 295. Р. 182-264.
3. Ivanchenko M.V. and Flach S. Anomalous conductivity: impact of nonlinearity and disorder // arXiv:1009.3447.
4. Flach S., Ivanchenko M.V. and Kanakov O.I. q-breathers in Fermi-Pasta-Ulam chains: Existence, local-
ization and stability // Phys. Rev. E. 2006. V. 73. Р. 036618.
5. Ivanchenko M.V. q-breathers and thermalization in acoustic chains with arbitrary nonlinearity index // Письма в ЖЭТФ. 2010. Т. 92. С. 405.
6. Lyapunov M.A. The General Problem of Stability of Motion. London: Taylor & Francis, 1992.
7. Гергель В.П., Стронгин Р.Г. Основы параллельных вычислений для многопроцессорных вычислительных систем // Учебное пособие. Нижний Новгород: Изд-во ННГУ им. Н.И. Лобачевского, 2003.
ANALYSIS OF EXACT PERIODIC SOLUTIONS IN LARGE OSCILLATORY CHAINS USING
PARALLEL COMPUTING
M.V. Ivanchenko
A parallel computing algorithm has been developed for numerical calculation and analysis of linear stability of exact periodic solutions in chains of nonlinear oscillators. The achieved acceleration has, for the first time, allowed such solutions to be obtained in large chains of more than one thousand elements, and also to show their localization in the space of system normal modes, to study their stability, and to compare the results with theoretical predictions. The practical significance of the results obtained is discussed.
Keywords: nonlinear dynamics, parallel computing, ensemble dynamics.
