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

УДК 519.622, 004.021, 004.942 Процедура построения симплектических численных схем для решения гамильтоновых систем уравнений

Б. Батгэрэл*^, Э. Г. Никонов*, И. В. Пузынин*

* Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, г. Дубна, Московская область, Россия, 141980 ^ Монгольский государственный университет науки и технологии Улан-Батор, Монголия

Применяемые в настоящее время схемы численного решения систем уравнений динамики многих частиц могут иметь ограничения по шагу и интервалу интегрирования, поскольку при их увеличении численные схемы становятся неустойчивыми и не сохраняют существующие интегралы движения.

В результате при моделировании поведения системы частиц на достаточно больших интервалах времени приходится уменьшать шаг интегрирования, что приводит к существенному увеличению объёма вычислений. В работе предложена новая процедура построения симплектических численных схем для решения гамильтоновых систем уравнений. Предложен подход к симметризации полученных симплектических разностных схем. Численные схемы, построенные при помощи предложенной в работе процедуры, сохраняют энергию системы на большом интервале численного интегрирования при относительно больших величинах шага интегрирования по сравнению с методом Верле, который обычно используется при решении уравнений движения молекулярной динамики. Приведены результаты численных экспериментов, демонстрирующие основные преимущества полученных в работе симметричных симплектических численных схем третьего порядка точности по шагу интегрирования гамильтоновых систем уравнений по сравнению с численными схемами метода Верле второго порядка.

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

1. Введение

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

Существующие пакеты программ молекулярной динамики для численного интегрирования уравнений движения используют схемы метода Верле. В них интегрирование ведётся с достаточно малым шагом по времени и контролем сохранения гамильтониана (энергии системы), поскольку уже при значении приведённой величины шага порядка 0.1 и количестве шагов порядка 103 накопленная вычислительная ошибка приводит к неустойчивости численной схемы и потере сохранения гамильтониана.

Статья поступила в редакцию 20 февраля 2016 г. Работа выполнена при поддержке гранта РФФИ 15-01-06055а.

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

— усовершенствования существующих пакетов:

— векторизация, как например, в пакете БЬ-РОЬУ [2];

— распараллеливание и ускорение межпроцессорных обменов на многопроцессорных системах;

— разработка спецпроцессоров [3];

— разработка эффективных новых и оптимизация существующих методов численного интегрирования уравнений движения.

Авторами предлагается новый подход к построению численных схем для решения уравнений движения молекулярной динамики. Будет показано, что построенные численные схемы сохраняют в приближенном решении основные свойства точного решения и обладают более высоким порядком аппроксимации, чем схемы метода Верле, а также являются устойчивыми по отношению к шагу аппроксимации.

Предложенный подход состоит в следующем:

— использование гамильтоновой формулировки уравнений движения молекулярной динамики;

— разложение точного решения в ряд Тейлора [4]. При этом возможно применение аппарата компьютерного алгебры для получения аналитических выражений для производных ряда Тейлора;

— использование для вывода численных схем аппарата производящих функций для сохранения геометрических свойств точного решения [1].

В данной работе решаются следующие задачи:

— построить симметричные симплектические численные схемы интегрирования гамильтоновых систем уравнений более высокого порядка точности, чем схемы метода Верле, что позволит увеличить шаг при фиксированном интервале интегрирования и, соответственно, уменьшить общий объем вычислений;

— выполнить численные эксперименты для сравнительного анализа свойств схем метода Верле и построенных численных схем;

— оценить перспективы применения построенных численных схем в существующих пакетах, например в ЬРМБ [5].

2. Постановка задачи

Движение системы N материальных точек в поле с потенциалом V(д), где q = — координата частицы, d = 3М — размерность пространства

координат, может быть описано с использованием системы гамильтоновых уравнений

^ (1)

дН(р, д) ^

дН(р, д) др

(2)

с начальными условиями

р(0) = р0, д(0) = д0

(3)

Здесь р = (р1,... ,ра)Т — импульс частицы, Н(р, д) — гамильтониан системы:

Н(р, д) = 1 рТМ(д) р + V(д),

(4)

М(д) — симметричная и положительно определённая матрица масс.

Гамильтонова система уравнений движения эквивалентна уравнениям, полученным в рамках ньютоновского формализма, если силы, действующие на материальные точки, представить в виде Г(д) = УУ(д) и стандартным образом

перейти от дифференциального уравнения Ньютона второго порядка заменой переменных к системе уравнений первого порядка

р = -f (Ч), (5)

<1 = р.

3. Численные методы

3.1. Геометрические методы

Геометрическим называют такой метод, который сохраняет некоторые геометрические свойства точного решения системы (1)—(3) [1].

Можно выделить следующие геометрические свойства решений:

- отображение ^ : (р^0),д(£0)) ^ (р^),д(1)), реализующее решение гамильто-новой системы уравнений (1)—(3), является симплектическим;

- решение обратимо во времени;

- решение сохраняет значение гамильтониана (4) для любого момента времени.

3.2. Симплектические методы

Дифференцируемое отображение

(р : и ^ (и е

называется симплектическим [1], если якобиан ^'(р, ч) удовлетворяет тождеству

Ч>'(Р, Ч)Т J Ч>'(Р, Ч) = з, (6)

где

т =( 0а

3 = V -I*

0^, - нулевая и единичная матрица размерности <1.

3.3. Симплектические численные методы

На дискретном множестве {Ък : к = 0,1,...; — Ьк = И} одношаговый метод решения системы (1)-(3), при постоянном шаге по времени И, можно представить в виде

где Ф^(р, ч) - преобразование приближенного решения (р% чfc)

при Ь — Ьк в приближенное решение (р^1, Ч^1) при £ =

Одношаговый метод (7) называется симплектическим, если преобразование (7), реализующее приближенное решение гамильтоновой системы (1)-(2), является симплектическим. В качестве примеров симплектических численных методов 1-го и 2-го порядка аппроксимации можно привести явный и неявный методы Эйлера и метод Верле. Метод Верле:

И2

Ч^+1 = <к + Ирк — у ¥ (чfc), (8)

р"+1 = р" — И ¥(1*) + ¥(Ч"+1)] . (9)

(р^+1, ч^+1) = ф^(р^, ч^), (7)

Алгоритм вычисления метода Верле:

р*+1/2 = р* - -Г(д*), (10)

д^1 = д^ + -р^1/2, (11)

р^+1 = р^+1/2 - - £ (д^+1). (12)

Как уже упоминалось, метод Верле является основным методом численного интегрирования уравнений движения в молекулярной динамике.

Если переписать систему уравнений (5) в виде обыкновенного дифференциального уравнения

д= - (д), (13)

то метод Верле может быть записан в форме стандартной разностной схемы второго порядка

д -2—2 +д = -г(д*), к = 0,1,.... (14)

Большинство стандартных численных методов изначально не являются сим-плектическими, в частности методы Рунге-Кутты. Для того, чтобы метод стал симплектическим, необходимо проделать ряд нетривиальных шагов для модификации метода в сторону симплектичности.

3.4. Симметричные численные методы

Одношаговый метод (7) называется симметричным, если удовлетворяет условию

Фь = Ф-£. (15)

В качестве примера симметричного симплектического метода можно привести метод средних точек, который для гамильтоновых систем вида (1), (2) может быть

записан следующим образом

р<=+. = р<= _ (р" +р"+1, ), (1в)

= д» + (р* +2рЬ+1 , ^^) . (17)

4. Симплектические преобразования и производящие

функции

4.1. Канонические преобразования

В гамильтоновой механике каноническое преобразование - это преобразование канонических переменных и гамильтониана

р :(р, д)^ (Р, д) (18)

не меняющее общий вид уравнений Гамильтона [в].

Канонические преобразования взаимнооднозначно определяются производящей функцией 5(р, д), которая является решением уравнения Гамилтона-Якоби и полный дифференциал которой равен

dS = Рт^д - рт^д.

(19)

4.2. Производящие функции

Производящая функция может быть выражена через любую пару из четырёх переменных р, q, Р, д. Возможны четыре варианта выбора пар переменных. Получаемые при этом функции принято называть производящими функциями 1-го, 2-го, 3-го или 4-го типа соответственно.

Производящие функции Производные

с = 5^, д, ¿) р = _ двг р = двг р дя Р дц

с = &(,, р, ¿) р = 052 д = р = ая д = ар

с = 5з (р, д, г) , = _ Мз Р = _ М3 , = ар Р = эц

с = 54(р, р, г) , = Ма. д = _ ад , = ая д = ар

4.3. Связь симплектических преобразований и производящих

функций

В соответствии с известной теоремой Пуанкарэ [1] для любого симплектиче-ского преобразования ^ существует производящая функция 5(,, д), и наоборот, для конкретной производящей функции 5(,, д) существует симплектическое преобразование которое может быть реконструировано при помощи следующих формул

яс яс

р = -л,^ "»• Р = ад ^ <20)

5. Процедура построения симплектических разностных схем на основе производящей функции первого типа

Для получения процедуры построения симплектических разностных схем была выбрана производящая функция первого типа 5 = 5*1 (,, д, ¿). При этом могут быть получены два варианта симплектических численных схем. Первый вариант получается при разложении канонической переменной , при £ = ¿&+1 = + Н в ряд Тейлора по переменным , и р при £ = , так называемом разложении в ряд Тейлора "вперёд". Второй вариант получается при разложении канонической переменной , при £ = в ряд Тейлора по переменным , и р при £ = так называемом разложении в ряд Тейлора «назад». При этом получаются соответственно явная и неявная схема до 3-го порядка включительно. Явно-неявная схема получается при т = 3 только в случае разложения в ряд Тейлора «назад». Причём все эти схемы изначально симплектические. В случае же т > 3 при использовании разложения в ряд Тейлора и «вперёд» и «назад» получаются явно-неявные симплектические разностные схемы.

В дальнейшем при построении численных схем использованы обозначения:

¿ = ^, р = р', , = q^

* = ^+1 = и + н, Р = р'+1, д =

Для производящей функции используются следующие обозначения: вт,1 — в случае разложения в ряд Тейлора «вперёд» и Ят,2 — в случае разложения в ряд Тейлора «назад». Здесь т — порядок аппроксимации точного решения.

5.1. Схема с разложением в ряд Тейлора «вперёд»

1. Представить в виде разложения в ряд Тейлора в точке с точностью

до 0(Нт+1).

2. Заменить pfc в разложении разностной производной от ч, кроме второго члена разложения. Получим уравнение относительно ч^1.

3. Решить полученное в предыдущем пункте, вообще говоря, нелинейное, для схем выше 3-го порядка аппроксимации, уравнение относительно ч^1. В случае нелинейного уравнения можно использовать метод Ньютона.

4. Выразить pfc через qfc и ч^1 с помощью полученного в первом пункте разложения в ряд Тейлора.

5. Найти производящую функцию Бт,1 = 5'(qfc, ч^1 ) путём интегрирования

и

полученного в предыдущем пункте выражения по .

6. Найти р^1 с помощью равенства

рк+1 = дБт,1

р = д ч^+1 ()

путём дифференцирования.

5.2. Схема с разложением в ряд Тейлора «назад»

1. Представить qfc в виде разложения в ряд Тейлора в точке ¿&+1 с точностью до 0(ЬГ+1).

2. Заменить р^1 в разложении разностной производной от ч, кроме второго члена разложения.

3. Выразить р^1 через qfc и ч^1 с помощью полученного выше разложения в ряд Тейлора.

4. Найти производящую функцию Бт,2 = 5'(qfc, Ч^1) путём интегрирования полученного в предыдущем пункте выражения по ч^1 .

5. Найти pfc помощью равенства

^ = — ^ (22)

путём дифференцирования.

6. Решить, вообще говоря, нелинейное относительно ч^1 уравнение (22). Везде в дальнейшем используется метод Ньютона.

7. Получить р^1 при помощи полученного во втором пункте выражения, ис-

и

пользуя известное значение и полученное в предыдущем пункте значение

5.3. Симплектическая разностная схема 2-го порядка

Для вычисления величины pfc используем разложение для < в ряд Тейлора «вперёд» до второго порядка

И2

Ч^1 = qfc + крк — у ¥ ) + 0(И3). В результате получим явную симплектическую разностную схему 2-го порядка

р* = ^^ + И ¥(Чк), (23)

^ = И ( Ч"+1И— Ч*) 2 — И [ V(<*) + V(ч^1) ] , (24)

(25)

Если для вычисления величины р'г+1 использовать разложение для q в ряд Тейлора «назад» также до второго порядка

Ъ2

qk = qfc+1 _ Ърк+1 _ ^{^^ + ^3),

то получим ту же самую схему (23), (25). Если же ввести обозначение

qk+1 _ qк __ к+1/2

Ъ

р

то получим известную схему Верле (10)—(12), которая является симплектической и симметричной.

5.4. Явная и неявная симплектические разностные схемы 3-го

порядка

Для построения явных и неявных схем третьего порядка, с использованием производящей функции первого типа, для нахождения величин рк и рк+1, используем разложения точного решения для канонической переменной q в ряд Тейлора для явной схемы «вперёд», для неявной - «назад».

Явная схема.

Ъ2 Ъ3

qk+1 = qk + Ърк _ г(qk) _ ^УГ(qk)

2

6

р

+ о(ъ4)

Если заменить рк в последнем члене ряда разностным отношением

Чк+1 _ г*к

q

р

+ О(К),

то получим явную схему 3-го порядка к _ qk+1 _ qk Ъ

^3,1 — —

р

Ч 2 _ Ъ[ 2у (qk)+у ^к+1) ] _ г )

Ъ

+ 2 г ) +у у г (qk)

^ qk+1 _ qk

Ъ

ъ2 с/ к\ qk+1 _ qk 6

Ъ

р

"+1 - qk+1 qk _ Ъ [ Г) + 2Г ^к+1) ] .

Ъ

Введя обозначение

0А

Ъ2

I + - УГ ^к)

-1

получим формулы для вычисления р и q на следующем временном шаге

(26)

(27)

(28)

к+1/2 — р к

Ъ.

р '- — рк _ 2Г(qk), р к+1/2 — @к рк+1/2,

2

= + Нр -+1/2,

р'+1 = р-+1/2 _ Н [ f (,'+1) ] . 6

Неявная схема.

Для получения неявной схемы представим значение для переменной , на й-ом шаге в виде разложения в ряд Тейлора «назад»

Н2 . Н3

= _ Нр'+1 _ — f(,'+1) + — Vf (,'+1) ■ р'+1 + 0(Н4)

,-+1

2 ^ ' 6

Представим р'+1 в последнем члене разложения в виде разностной производной

Ч-+1 _

р , - = _

Н

+ О(Н).

В результате получим следующее выражение для р'+1:

р'+1 = ,"+1Н_ _ Нf (,"+1) + Н2 Vf (,"+1) ■ ,"+1Н_ ^ . (29)

После интегрирования по получаем производящую функцию £3,2

7 / ' + 1 — \ 2 7 7 2 '+1 —

5з,2 = 2 (_ Н [ ^(,к) + 2^(,к+1) ] + Нтf(,к+1) ■ (30)

и дифференцированием по с учётом знака находим выражение для р-:

р- = ,-+1н_ + Н [ 2f ) + f (,-+1) ] . (31)

Далее аналогично явной схеме получаем формулы для неявной схемы третьего порядка

р-+1/2 = р- _ Нf), (32)

,-+1 = + Нр-+1/2, (33)

Р(,-+1) = ,-+1, (34)

Д-+1 = Н [ f (,-+1) _ f (,-) ] , (35)

Н2

0-+1 = I + — Vf (,-+1), (36)

6

р-+1/2 = 0-+1(р-+1/2 _ д-+1), (37)

р-+1 = р-+1/2 _ Нf(,-+1), (38)

где функция Е(,-+1) определяется по формуле

Н2

Г(,-+1) = ,-+1 + _ [ f (,-+1) _ f (,-) ] . (39)

6. Симметричные симплектические разностные схемы

Для построения симметричных симплектических разностных схем рассмотрим одно параметрическое семейство производящих функций следующего вида.

Ят(а) — аЯт,1 + (1 _а)Ят,2. (40)

Здесь а : 0 < а < 1, Ят,1 и Ят,2 - производящие функции, полученные при помощи процедур (21) и (22) соответственно. Каждая из этих производящих функций Бт(а) порождает симплектическую разностную схему т-го порядка. Если а — 1/2, то мы получаем процедуры построения симметричных симплектических разностных схем.

Получим процедуры построения симметричных симплектических разностных схем 3-го и 4-го порядка, поскольку схемы 3-го и 4-го порядка представляют наибольший практический интерес при проведении численных расчётов. Для простоты представления схема 4-го порядка приведена для одномерного случая. Схемы более высокого порядка аппроксимации имеет смысл выводить с использованием средств компьютерной алгебры.

6.1. Симметричная симплектическая схема для т — 3

Ъ / qk+1 _ qk \ 2 Ъ

^з,з — 2(_ Ъ (ч-) + V] +

+ [ Г_ Г) ] ■ ^ qk , (41) р к — q-+1p- + 12 [ (qk) + Г(qk+1) ] + УГ(ч-) ■ q-+1Ъ-^-, (42) рк+1 — q-+1p- _ 12 [ г(qk) + (qk+1) ] + УГ(qk+1) ■ . (43)

6.2. Симметричная симплектическая схема для т — 4

Для построения симметричной симплектической численной схемы 4-го порядка сначала проведём вывод схемы «вперёд». Сначала представим Як+1 в виде разложения в ряд Тейлора «вперёд»

Ъ2 ъз ъ4

Як+1 — Як + Ърк _ — ДЯк) _ у //(Як)рк _ 24 [/"(1к)рк _ //(ЯкЖЯк)] + о(ъЪ). После замены рк в члене разложения с Ъ3

Як+1 _ Як . Ъ /п/ь2\

Рк — -Ъ- + ^Д 1к) + 0(П )

и замены квадрата рк в члене разложения с Ъ4

к+1 _ к

к

Ъ

+ О(Ъ)

получим следующее выражение для рк, зависящее только от о- и Як+1'-

+

Як+1 _ Як + Ъ )+ Ъ2 ) Рк —-Ъ-+ ^Д Чк) + у / (Як)

Як+1 _ Як Ъ

+ 2/(<?к)

Нз + 24

/"(9-)(2_ /'(9-)/( ,-)

Далее интегрированием получаем производящую функцию

Н /д-+1 _ дЛ _ Н 2 V Н У 4

т2 . - _ /7, \ Т^3 — _ /7, \ 2 Ъ3

54,1 = £( _4[3У (9- ) + У (,-+!)]_

\ ^-+1 _ Н3 гЧ \ (^-+1 _ Н3 Г ¿2/ \ ^^ М _ ^Д9-) ( -Н-) _ 24 ^ (Н-Н-У _ 48 ^ (9-+1) _ ^ () ^

4' V Н У 24' V Н У 48 и дифференцированием выражение для -+1

Р-+1 = ^-+1Н ^ _ Н [ /(9-) + Я9-+0 ] _

_ у^/Ч9-) (^-+1Н ^ + 724/Ч9-+1)^(9-+1).

Аналогично, в соответствии с процедурой (22), строится схема «назад» 4-го порядка. Сначала представим д- в виде разложения в ряд Тейлора «назад»

н2 Н3

= 9-+1 _ НР-+1 _ Ну/(9-+1) + у/'(<?-+!)р-+1 _

74

_ [ /''(9-+1)р-+1 _ /'(9-+1)/(9-+1) ] + ^(Н5).

3

Затем проведём замену р-+1 в члене разложения с Н

-+1 - Н 2 №+1 =-у--) + °(Н ^

и замену квадрата -+1 в члене разложения с Н4

-+1 -Р-+1 = —у--+ О(Н).

В результате получим выражение для р-+1, зависящее только от д- и ^-+1

3-+1 _ Н „ )+ Н2 и ,-+1 ,- ,+ р-+1 = —у--9-+1) + (9-+1)' —^— ' +

Н3

+ 24

/ 9-+1 _ О- \

V Н / /"(?-+!)(2 +/'(,-+1)/(9-+1)

Далее интегрированием получаем соответствующую производящую функцию 54,2 = Н (^р-)2 _ 4 [ у (?-) + ЗУ (,-+1) ] +

/ 0-+1 _ д- \ V Н У

Н2 -+1 - Н3 -+1 -

+т/( ^Ч ""н- ) _ ^Ч ""~)

н3

+ Н8 [ /2(9-+1)_ /2(9-)]

и дифференцированием выражение для р- :

Як+1 _ Чк , Ъ л N1 Ъ2 ГЦ \(Як+1 _ Чк\ . ъ3

рк

ъ + Ъ [ Я) + /(Як+1) ] _ ^/Ч9к+1Ъ ^ + т;^/'(о-)/( Далее с помощью формулы (40) при а — 0.5

$4,3 — 1 (V + ^4,2) получаем симметричную производящую функцию 54,3

,4,3 — Ъ ( )- ^ (. к) + , (.-+,)]+Ъ2 [ /(а+1) _/(.-)]

_ Ъ [ ^ <®)+^ (®+1) ]+Ъ- [ /(«+1) _ /(«) ] ( "-'^р- ) _ [ /'(®+1)+/'(») ] (2+48 [ л «-+■) _ /2( ) ] ■

|[/'(«+1)+л«)](2+> -

В результате после соответствующей процедуры дифференцирования получим симметричную симплектическую схему 4-го порядка:

р-—+|[„ы+/(. к+1)]_ |4 [2,,-) , +

+ Ъ [ 3 /( ?к) + /(^+1) ] _ Ъ4 [ 2 АЯк) _ /'(?к+1) ] ()

+ Ъ8/''() (^р- )2 _ Ъ4/'(®)/( ),

рк+1 — ^^_Ъ [ Кяк) + з/(,к+1) ]_Ъ4 [ А9к) _ 2/4^+1) ] (^^)

7. Численные эксперименты

Численные эксперименты выполнены для сравнительного анализа численной схемы метода Верле и численных схем более высокого порядка аппроксимации при т — 3, полученных с использованием разработанной авторами процедуры, на примере задачи Кеплера [1]. В ней движение двух материальных точек описывается гамильтонианом следующего вида

н(р,я) — 2(р1 + р) _ ■ (44)

2 V <?1 + ъ

Численные расчёты проводились при следующих параметрах:

1) Ъ — 0,01, Т — 1000;

2) Ъ — 0,05, Т — 1000;

3) Ъ — 0, 10, Т — 1000;

4) Ъ — 0, 20, Т — 5000.

Здесь Ъ - шаг по времени, Т - граница интервала времени [0,Т] эволюции системы.

Для каждого фиксированного шага по времени Ъ сравнивались результаты расчётов, полученные с использованием метода Верле, и результаты, полученные

с использованием построенных в работе явной, неявной и симметричной численной схемы 3-го порядка. В качестве «точного» решения использовались результаты расчётов по методу Верле с шагом по времени, равным Л./10. Результаты моделирования представлены на рис. 1-14. Задача Кеплера. Явная схема.

The Kepler 2 body problem Time = 1000

The Kepler 2 body problem Hamiltonian

^flhtfafflMMf^ek^fw

■i ¡1 s ji A' 4

ii

-m3A time step= 0 01 ■1 Verlet time step= 0 01 — Verlet time step= 0 001

У"?

!

0 4 \ /

I

l' 'l Q '' '! П a " f"

а) Фазовые траектории б) Зависимость гамильтониана

от времени

Рис. 1. Явная схема, схема Верле при к = 0.01 и схема Верле при к = 0.001

The Kepler 2 body problem Hamiltonian

а) Фазовые траектории Рис. 2. Явная схема, схема Верле при к = 0.05 и схема Верле при к = 0.005

б) Зависимость гамильтониана от времени

The Kepler 2 body problem Time = 1000

The Kepler 2 body problem Hamiltonian

-2 -1 5

'1 Ü Ii i" ü ! i; 1 i! J ü ii II ~-i 1 -- i: с ! i < • »

8 ! — H — m3A time step= 0 1 j[ - i ^^ i Verlet time step= 0 1 | ! -Verlet time step= 0 01 [

•Iii* A 4 Тф

• > -; ' '11 у 1 ¡is«.-, ..^»iM':^^, r.if

а) Фазовые траектории Рис. 3. Явная схема, схема Верле при к = 0.1 и схема Верле при к = 0.01

б) Зависимость гамильтониана от времени

200

400 600

800

000

qi

I

200

400 600

800

000

qi

0 qi

200

400 600

800

000

The Kepler 2 body problem Time = 600

The Kepler 2 body problem Hamiltonian

-m3A time step= 0 2

- • Verlet time step= 0 2

- Verlet time step= 0 02

а)

Рис

Фазовые траектории б) Зависимость гамильтониана

от времени

. 4. Явная схема, схема Верле при к = 0.2 и схема Верле при к = 0.02

-4

-4

q1

Как показывают численные эксперименты, при величинах шага по времени до Ъ — 0.1 включительно схема Верле и явная схема третьего порядка очень близки и по расчётным фазовым траекториям и по значениям гамильтониана на всем интервале интегрирования (рис. 1-3). Начиная с величины шага Ъ — 0.2 схема Верле оказалась более устойчивой, хотя и менее точной по сравнению с явной схемой третьего порядка. Кроме того, начиная примерно с 500-го шага она становится неустойчивой (рис. 4) по шагу интегрирования. Задача Кеплера. Неявная схема.

The Kepler 2 body problem Time = 1000

The Kepler 2 body problem Hamiltonian

а) Фазовые траектории

б) Зависимость гамильтониана от времени

Рис. 5. Неявная схема, схема Верле при к = 0.01 и схема Верле при к = 0.001

The Kepler 2 body problem Time = 1000

The Kepler 2 body problem Hamiltonian

-0 492 -0 494

-0 498 J*

-I*

-m3B timestep=0 05 ■' Verlet time step= 0 05 -Verlet time step= 0 005

Ж

X

¥

»10

. i 'i?

а) Фазовые траектории б) Зависимость гамильтониана

от времени

Рис. 6. Неявная схема, схема Верле при к = 0.05 и схема Верле при к = 0.005

qi

-0 5 0

qi

200

400 600

800

000

The Kepler 2 body problem Time = 1000

а) Фазовые траектории

The Kepler 2 body problem H

— m3B time step= I Verlet time step= Verlet time step=

'l.l'lp 'l .1 " n T

i|T

rh

1 " I 4 Ï 4 i' I « . □

I11"' ÎS°9Br, BD<o!i '< 1 IV°o

USM4«<

б) Зависимость гамильтониана

от времени

200

400

600

800

000

qi

Рис. 7. Неявная схема, схема Верле при h = 0.1 и схема Верле при h = 0.01

The Kepler 2 body problem Time = 320

-m3B time step= 0 2 -■Verlet time step= 0 2 -Verlet time step= 0 02

а) Фазовые траектории Рис. 8. Неявная схема, схема Верле при к = 0.2 и схема Верле при к = 0.02

б) Зависимость гамильтониана от времени

6 -Ъ -4 -3 -2 -1 0 1 2 3

qi

Аналогично численным экспериментам с явной схемой третьего порядка, в случае неявной схемы, при величинах шага по времени до Н = 0.1 включительно, схема Верле и неявная схема третьего порядка очень близки и по расчётным фазовым траекториям и по значениям гамильтониана на всем интервале интегрирования (рис. 5-7). Начиная с величины шага Н = 0.2 схема Верле снова оказывается более устойчивой по сравнению с неявной схемой третьего порядка, однако примерно с 250-го шага становится неустойчивой по шагу интегрирования (рис. 8).

Задача Кеплера. Симметричная схема.

The Kepler 2 body problem Time = 1000

а) Фазовые траектории

б) Зависимость гамильтониана от времени

Рис. 9. Симметричная схема, схема Верле при h

h = 0.001

0.01 и схема Верле при

200

400

600

800

300

qi

The Kepler 2 body problem Hamiltonian

а) Фазовые траектории

б) Зависимость гамильтониана от времени

Рис. 10. Симметричная схема, схема Верле при к

к = 0.005

0.05 и схема Верле при

The Kepler 2 body problem Time = 1000

The Kepler 2 body problem Hamiltonian

■i 11

i Hl 1 i ■■ r I i ■ i 4

■ — В — m3C time step 1 i ^^ i Verlet time step -Verlet time step 0 1 0 1 0 01 i" • 6

i 1 i , 4

в In ------ "Iii

а) Фазовые траектории

б) Зависимость гамильтониана от времени

Рис. 11. Симметричная схема, схема Верле при к = 0.1 и схема Верле при

к = 0.01

The Kepler 2 body problem Time = 1000

lisifai

ix MVS ........... 11 - ' ^

The Kepler 2 body problem Hamiltonian

— m3C time step= 0 2 i

— • Verlet time step= 0 2 ;

— Verlet time step= 0 02

а) Фазовые траектории.

б) Зависимость гамильтониана от времени

Рис. 12. Симметричная схема, схема Верле при к = 0.2 и схема Верле при

к = 0.02

Совсем иная картина открывается при сравнении результатов расчётов с использованием симметричной симплектической численной схемы третьего порядка. Вплоть до величины шага численного интегрирования по времени Н = 0.2 включительно симметричная симплектическая схема третьего порядка существенно более точна по расчётным фазовым траекториям и по отклонениям значений

I

200 400 600 800

000

qi

2UU

4UU 600

800

UUU

-4

-4

0 q1

гамильтониана от начального значения на всем интервале интегрирования (рис. 9-12).

Задача Кеплера. Симметричная схема. Т = 5000.

The Kepler 2 body problem: Time = 5000

The Kepler 2 body problem: Hamiltonian

---m3C time step= 0.2

-----Verlet time step= 0.2

— Verlet time step= 0.02

а) Фазовые траектории

2000 3000 4000

time

б) Зависимость гамильтониана от времени

Рис. 13. Симметричная схема, схема Верле при к = 0.2 и схема Верле при

к = 0.02

14

-4

0

2

4

0

1000

5000

qi

Number of iteration: e = 10 10

Рис. 14. Количество итераций метода Ньютона на каждом шаге по времени

при £ = 10-10

И, наконец, в численных расчётах, которые проводились на большем интервале эволюции системы Т = 5000, симметричная симплектическая численная схема третьего порядка продемонстрировала существенно большую точность и устойчивость при величине шага Н = 0.2. При этом схема Верле становится неустойчивой уже начиная с интервалов порядка Т = 4000 (рис. 13). При это вычислительные затраты в методе Ньютона при вычислениях для неявной части симметричной схемы не превышают 3 итераций при точности вычислений е = 10-10 (рис. 14).

S. Заключение

1. Предложен подход к построению симметричных симплектических численных схем интегрирования уравнений движения метода молекулярной динамики в гамильтоновой формулировке.

2. Описан алгоритм получения симметричных симплектических разностных схем заданного порядка аппроксимации с использованием аппарата производящих функций.

3. Приведены симметричные симплектические разностные схемы до 4-го порядка аппроксимации включительно. Схемы более высокого порядка аппроксимации требуют более сложного вывода и в дальнейшем могут быть получены с использованием аналитических вычислений.

4. Проведены численные эксперименты, показавшие, что полученные на основе разработанного подхода разностные схемы третьего порядка аппроксимации сохраняют устойчивость для достаточно больших значений шага и интервала интегрирования в отличие от метода Верле. Они с более высокой точностью сохраняют гамильтониан системы на всем интервале интегрирования.

Литература

1. Hairer E., Lubich C., Wanner G. Geometric Numerical Integration. — Berlin, Heidelberg: Springer, 2006. — P. 644.

2. Методы молекулярной динамики для моделирования физических и биологических процессов / Х. Т. Холмуродов, М. В. Алтайский, И. В. Пузынин и др. // Физика элементарных частиц и атомного ядра. — 2003. — Т. 34, № 2. — С. 472515.

3. MDGRAPE-4: a Special-Purpose Computer System for Molecular Dynamics Simulations I I. Ohmura, G. Morimoto, Y. Ohno et al. // Philos Trans A Math Phys Eng Sci. — 2014. — Vol. 372, No 2021. — Pp. 1-16.

4. Akishin P. G., Puzynin I. V., Vinitsky S. I. A Hybrid Numerical Method for Analysis of Dynamics of the Classical Hamiltonian Systems // Comp. Math. Applic. — 1997. — Vol. 34, No 2-4. — Pp. 45-73.

5. Las Palmeras Molecular Dynamics: Flexible and modular molecular dynamics / S. Davis, C. Loyola, F. Gonzalez, J. Peralta // Computer Physics Communications. — 2010. — Vol. 181, No 12. — Pp. 2126-2139.

6. Гантмахер Ф. Р. Лекции по аналитической механике. — М.: Физматлит, 2001. — С. 262.

UDC 519.622, 004.021, 004.942

A Procedure for Constructing Simplectic Numerical Schemes for Solving of Hamiltonian Systems of Equations

B. Batgerel*t, E. G. Nikonov*, I. V. Puzyninf

* Laboratory of Information Technologies Joint Institute for Nuclear Research 6, Joliot-Curie, Dubna, Moscow region, Russia, 141980 ^ The Mongolian University of Science and Technology 8th khoroo, Baga toiruu 34, Sukhbaatar district Ulaanbaatar, Mongolia, 14191

Numerical schemes which is used for solving of many-particle dynamics systems of equations can have restrictions on a step and an interval of integration because if its increase the numerical schemes became unstable and don't conserve existing integrals of motion. As a result when we simulate many-particle system behavior on the sufficiently large time interval we should decrease an integration step which leads to considerable increasing of computation quantity. In this paper a new procedure for constructing simplectic numerical schemes for solving of Hamiltonian systems of equations is proposed. A method for symmetrization of received simplectics numerical schemes is proposed too. Constructed by proposed in the

paper procedure numerical schemes conserve energy of a system on the large interval of numerical integration for relatively large integration step in comparison with Verlet method which is usually used for solving of equations of motion in molecular dynamics. Results of numerical experiments are given in the paper. These results show main advantages of received symmetric simplectic numerical schemes of third order of accuracy for the integration step for the Hamiltonian systems of equations in comparison with numerical schemes of Verlet method of second order of accuracy.

Key words and phrases: Hamiltonian systems of equations, simplectic difference schemes, generating functions, molecular dynamics.

References

1. E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration, Springer, Berlin, Heidelberg, 2006.

2. H. T. Holmurodov, M. V. Altajskij, I. V. Puzynin, et al., Molecular Dynamics Methods for Simulation of Physical and Biological Processes, Phusics of Elemrntary Particles and Atomic Nuclei 34 (2) (2003) 472-515, in Russian.

3. I. Ohmura, G. Morimoto, Y. Ohno, A. Hasegawa, M. Taiji, MDGRAPE-4: a SpecialPurpose Computer System for Molecular Dynamics Simulations, Philos Trans a Math Phys Eng Sci. 372 (2021) (2014) 1-16.

4. P. G. Akishin, I. V. Puzynin, S. I. Vinitsky, A Hybrid Numerical Method for Analysis of Dynamics of the Classical Hamiltonian Systems, Comp. Math. Applic. 34 (2-4) (1997) 45-73.

5. S. Davis, C. Loyola, F. Gonzalez, J. Peralta, Las palmeras molecular dynamics: Flexible and modular molecular dynamics, Computer Physics Communications 181 (12) (2010) 2126-2139.

6. F. R. Gantmakher, Lections on Analytical Mechanics, PhysMathLit, Moscow, 2001, in Russian.