Программирование, 2022, № 2, стр. 47-52

1. ВВЕДЕНИЕ/МОТИВАЦИЯ

Данная статья – продолжение серии работ, касающихся геометрических интеграторов, применяемых для компьютерного моделирования механических систем. Изучение геометрических интеграторов – численных методов сохраняющих внутреннюю геометрическую структуру дифференциальных уравнений – в свою очередь, является важной частью более глобального проекта по “геометризации механики”.

Наиболее хорошо изученным примером этого подхода является использование симплектических и вариационных интеграторов для консервативных конечномерных механических систем. Динамика системы описывается в непрерывном случае либо Гамильтоновыми уравнениями, внутренняя структура которых задается симплектической формой, либо уравнениями Эйлера–Лагранжа, полученными применением вариационного принципа к некоторому функционалу (см. например, [1]). При дискретизации этих уравнений в Гамильтоновом случае сохраняется симплектическая форма (см. обзор таких методов в [2]). Тогда можно показать [3], что, так как динамика задается в точности функцией полной энергии системы, некоторая ее дискретная версия будет тоже сохраняться – разницу между этой дискретной версией и исходной непрерывной функцией можно оценить и тем самым контролировать сохранение энергии в численном счете. В Лагранжевой картине строится дискретный аналог вариационного принципа (ДВП, [4]), что позволяет не только выписать соответствующую разностную схему, но и подобрать шаг интегрирования для точного сохранения энергии.

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

Под обобщенной геометрией мы понимаем, в первую очередь, структуры Дирака, которые естественным образом возникают при описании диссипативных и взаимодействующих систем, а также систем со связями. Последним и посвящена эта статья. В предыдущей работе по этой теме [6] мы сформулировали ряд открытых вопросов, для решения которых может быть полезна компьютерная алгебра. Один из них, а именно построение разностной схемы по заданной структуре Дирака, оказался концептуально более глубоким, чем просто применение существующих алгоритмов символьных вычислений – мы сформулируем его снова здесь и предложим алгоритмический метод его решения.

2. СОХРАНЕНИЕ СТРУКТУР ДИРАКА

Как мы уже упомянули, структуры Дирака возникают в достаточно общих ситуациях при качественном анализе механических систем. Они определены, в частности, с этой мыслью в оригинальной работе [7], упрощенное изложение которой мы привели в [6]. Для самодостаточности этой статьи напомним, что речь по-прежнему идет о двойственности описания механических систем: динамику можно рассматривать в терминах векторных или ковекторных полей на фазовом пространстве. Структуры Дирака дают геометрическое описание связи между двумя картинами: в каждой точке обобщенного фазового пространства выполняется некоторое линейное соотношение между объектами, а глобальная структура удовлетворяет дифференциальному условию инволютивности.

Структура Дирака для систем со связями – вполне естественный пример, с точки зрения геометра, она является (деформированным) прообразом так называемого распределения связей. Ее применение для построения численного метода, сохраняющего связи, по-видимому, впервые подробно описано в [8], основываясь на вариационном подходе [4]. В [9] мы построили естественную модификацию этого метода, позволяющую улучшить точность интегрирования. Именно этот процесс построения разностной схемы мы и предложили автоматизировать в [6].

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

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

Определение 1. Будем называть разностную схему псевдогеометрическим интегратором порядка ${\mathbf{(k,m)}}$, если она определяет численный метод порядка $k$ по решению и сохраняет геометрическую структуру с точностью до порядка $m$ по шагу интегрирования.

Понятно, что интересные случаи соответствуют значениям порядка $m > k$. В этих обозначениях настоящий геометрический интегратор будет псевдогеометрическим порядка $(k,\infty )$.

Ниже мы объясним, как такие интеграторы строятся для структур Дирака, избегая, в частности, необходимости “угадывания” дискретизации. Построение основано на использовании обобщения семейства методов Рунге–Кутта.

3. АЛГОРИТМ

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

Алгоритм 1.

1. Для всех динамических переменных системы записать метод Рунге–Кутта с числом стадий $s$ с неопределенными коэффициентами.

2. Предположить, что в момент времени $n$ все динамические переменные выполняют условие принадлежности структуре Дирака, оценить в общем виде погрешность этого условия в момент времени $(n + 1)$ с учетом приближения, записанного в пункте 1).

3. В порядке возрастания степеней разложения погрешности по шагу интегрирования приравнять к нулю максимально возможное количество слагаемых – получить условия на коэффициенты из пункта 1).

Несмотря на кажущуюся простоту данного алгоритма, важно понимать ряд деталей.

Замечание 1. Методы Рунге–Кутта, используемые в пункте 1) могут и, в общем случае, должны быть неявными, то есть таблицы Батчера будут полными, а не треугольными. Это наблюдение следует из того факта, что для частного симплектического случая структуры Дирака соответствующий интегратор должен воспроизвести известный симплектический [10].

Замечание 2. Для разных переменных методы могут и, в общем случае, должны иметь разные коэффициенты. Именно это обобщение стандартных методов Рунге–Кутта позволяет (ценой введения дополнительных динамических переменных) достичь дополнительной свободы, помогающей сохранять геометрическую структуру. Однако не нужно думать, что коэффициенты для разных переменных не будут повторяться совсем. Разумное разделение снова “подскажет” структура Дирака: естественно, например, объединить в три группы, соответственно, координаты, скорости и импульсы системы.

Замечание 3. Используемое условие на принадлежность структуре Дирака может содержать производные некоторых функций, характеризующих систему, но не производные динамических переменных – мы рассмотрим эту ситуацию в примере ниже. Это значит, что условие, выражаемое в пункте 2, чисто алгебраическое. Погрешность же получается явным его разложением в ряды Тейлора вокруг точки в момент времени $n$, используя пункт 1 как приращение аргумента.

4. ИМПЛЕМЕНТАЦИЯ

Каждый пункт приведенного алгоритма – совершенно явное символьное вычисление, которое разумно производить средствами компьютерной алгебры. Единственный этап, который мы не описали, связан непосредственно с нахождением коэффициентов из условий пункта 3. Заметим, что по построению, эти условия полиномиальные, а некоторые даже линейные. Поэтому задача решается стандартными методами компьютерной алгебры. Заметим также, что если решение уравнений на коэффициенты неединственно, разумно выбрать наиболее явный метод, чтобы ускорить дальнейший численный счет.

Основная цель данного алгоритма – построение разностных схем высокого порядка. То есть, число стадий методов Рунге–Кутта интересно сделать достаточно большим. С ростом числа стадий количество неопределенных коэффициентов растет квадратично, и промежуточные приближения будут входить в условия нелинейно – это может вызвать технические трудности связанные с решением системы из пункта 3. Напомним, однако, что для каждой группы переменных результат будет представителем семейства методов Рунге–Кутта. Это значит, что даже если самое общее решение получить не удается, например, из-за нехватки вычислительных ресурсов, можно рассмотреть комбинации известных классических решений.

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

5. ПРИМЕР СИСТЕМ СО СВЯЗЯМИ

Достаточно характерный пример применения приведенного выше алгоритма – изучение механических систем со связями. Речь идет о механических системах, на координаты и скорости которых наложены некоторые ограничения. Уравнения, описывающие их динамику имеют вид:

где $L(q,\dot {q})$ – Лагранжиан системы, а ${{\alpha }^{a}}$ – генераторы идеала порожденного связями, ${{\lambda }_{a}}$ – соответствующие множители Лагранжа. Условия на принадлежность этому идеалу обычно заданы явно и имеют вид ${{\varphi }^{a}}(q,\dot {q}) = 0$. На практике такие системы решаются введением дополнительных переменных¹¹ ${v}: = \dot {q}$, что позволяет переписать (5.1) как систему первого порядка. Систему можно численно интегрировать классическими методами, если множители Лагранжа явно вычисляются по координатам и скоростям в каждый момент времени.

Подход с использованием структур Дирака подразумевает, в дополнение к ${v}$, введение переменных $p$, соответствующим импульсам системы. Они получаются при рассмотрении Гамильтонова описания системы, то есть после преобразования Лежандра. Отличие от классического подхода заключается в том, что переменные ${v}$ и $p$ теперь будут рассматриваться одновременно, а их соотношение через преобразование Лежандра, а также связи, наложенные на систему, будут учитываться отдельными условиями – эти условия в точности соответствуют определению структуры Дирака ${{\mathbb{D}}_{\Delta }}$ для тройки $(q,{v},p)$ (см. [6]). Опуская геометрические подробности, приведем явный вид описанных выше условий:

здесь $\Delta $ обозначает распределение связей, то есть, в каждой точке конфигурационного пространства – множество совместных со связями скоростей, а ${{\Delta }_{0}}$, как и в уравнении (5.1), соответствующий идеал. Таким образом, при условиях (5.3), уравнения (5.2) эквивалентны (5.1) в непрерывном случае.

Рассмотрим для примера, согласно алгоритму выше, семейство двустадийных методов Рунге–Кутта с шагом $h$:

где в правых частях разностной схемы стоят стандартные приближения правых частей дифференциальных уравнений: и аналогично для $\tilde {l}$, $\bar {l}$, с очевидными изменениями. Подчеркнем еще раз, что методы могут быть неявными, то есть изначально все коэффициенты ${{a}_{{ij}}}$, ${{\tilde {a}}_{{ij}}}$, ${{\bar {a}}_{{ij}}}$ не равны нулю. Вспомним также, что все переменные системы многомерные, но коэффициенты мы различаем лишь для трех групп: для координат, скоростей и импульсов системы.

Теперь подставим разложения (5.4) в условия (5.3) на принадлежность структуре Дирака, которые в данном случае распадаются на два: импульсы определены преобразованием Лежандра: $\psi (q,p,{v}): = \left( {p - \frac{{\partial L}}{{\partial {v}}}} \right)$ = 0, и скорости удовлетворяют связям: ${{\varphi }^{a}}(q,{v}) = 0$. Выпишем, как и предполагает пункт 2 алгоритма соответствующие ряды Тейлора для $\psi $ и ${{\varphi }^{a}}$, с учетом уравнений движения (5.1) и (5.2). Приравняв к нулю коэффициенты при разных степенях $h$, получим следующее утверждение.

Теорема 1. Разностная схема (5.4) определяет численный метод по крайней мере второго порядка, если выполнены следующие условия на коэффициенты:

Этот численный метод сохраняет структуру Дирака ${{\mathbb{D}}_{\Delta }}$ и, как следствие, преобразование Лежандра и связи системы, по крайней мере с точностью до третьего порядка, если к тому же выполнены следующие условия:

а также:

Иными словами, при выполнении условий выше численный метод (5.4) является псевдогеометрическим интегратором по крайней мере порядка (2, 3).

Замечание 4. Первый набор условий в теореме (5.5) не определяет никаких зависимостей между разными группами переменных – это суть классические условия для методов Рунге–Кутта. Остальные условия разделены на две группы, поскольку уравнения (5.6) получаются из сохранения преобразования Лежандра, а (5.7) из сохранения связей, но для произвольной структуры Дирака, конечно, нет никаких причин рассматривать их отдельно.

Замечание 5. Формулировка “по крайней мере” для порядка метода не случайна. В данном примере легко понять, какие условия независимы. В общем же случае, когда рассматриваются многостадийные методы Рунге–Кутта продолжение разложения Тейлора может быть громоздким с сильно нелинейными зависимостями. То есть определить алгоритмически, получаются ли новые условия или следствия старых, может быть технически трудно или невозможно. Точное решение этого вопроса касается напрямую символьных вычислений и выходит за рамки данной статьи. На практике же он может решаться принудительной остановкой разложения до достижения максимального количества слагаемых в пункте 3) алгоритма.

Замечание 6. Заметим также, что достаточно часто уравнения (5.5)–(5.7) допускают множество решений, однако разложение до следующей степени шага дискретизации приводит к несовместным условиям. Тогда разумно рассмотреть частичное сохранение геометрической структуры: в нашем примере это могут быть только преобразования Лежандра или только уточнение связей.

Замечание 7. В приведенной теореме мы рассматриваем наиболее общий случай структуры Дирака для систем со связями. То есть мы предполагаем, что нет никаких известных функциональных зависимостей между $L$ и ${{\varphi }^{a}}$, и не ожидаем никаких упрощений рядов Тейлора из-за вырождения функций. В конкретных примерах механических систем некоторые условия могут быть ослаблены, или полученный метод может оказаться более высокого порядка.

В качестве иллюстрации метода мы применили данную теорему для системы описывающей движение математического маятника на плоскости в поле силы тяжести, подробно описанной в [9]. Динамика системы описывается в конфигурационном пространстве $Q = {{\mathbb{R}}^{2}}$, $q = ({{q}_{x}},{{q}_{y}})$, она задается свободным Лагранжианом

и связью $\varphi = {{\left\| q \right\|}^{2}} - {{l}^{2}}$. Структура Дирака ${{\mathbb{D}}_{\Delta }}$ задается стандартным условием $p = \frac{{\partial L}}{{\partial {v}}}$ и идеалом, порожденным $\alpha : = d\varphi $.

Таблица 1 показывает сравнение двух разных псевдогеометрических методов порядка (2.3) (RKD-2), построенных по нашему алгоритму, с двухшаговым методом из [9] (Dirac-2).

Из таблицы видно, что за счет неявности метода новый псевдогеометрический интегратор проигрывает в производительности, однако существенно выигрывает в вопросе сохранения геометрической структуры. Заметим, что даже метод Dirac-2 превосходил классические методы сравнимого порядка, а наш подход позволил значительно улучшить его точность.

Конечно, движение математического маятника – это лишь “академический пример”, где динамика хорошо изучена. Но именно поэтому он позволяет проиллюстрировать особенности интегратора, не отвлекаясь на, возможно, сложное нелинейное поведение системы.

6. ВМЕСТО ЗАКЛЮЧЕНИЯ

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

В рамках универсального подхода он позволяет воспроизвести результат [8] – предложенный там метод в наших обозначениях будет псевдо-геометрическим порядка (1.2) для структуры Дирака. Кроме того, когда структура Дирака строится по симплектической структуре (случай описанный здесь, но без связей), метод позволяет восстановить симплектический метод описанный в [11] и все семейство симплектических методов Рунге–Кутта [10].

С точки зрения приложений в геометризации механики, наш подход обладает важным преимуществом: он применим как для Лагранжева формализма, так и для Гамильтонова описания. Первый важен, например, в контексте вариационной формулировки динамики на структурах Дирака [12], второе естественно возникает при построении разностных схем, сохраняющих структуры Пуассона [13]. А для конкретных задач компьютерного моделирования применение псевдогеометрических интеграторов позволяет более точно различать динамические режимы даже для потенциально хаотических систем [14].