О вычислении переменных разделения в уравнении Гамильтона--Якоби на компьютере Текст научной статьи по специальности «Математика»

УДК 517.9+531.01

О вычислении переменных разделения в уравнении Гамильтона—Якоби на компьютере

Ю. А. Григорьев, А. В. Цыганов

Санкт-Петербургский государственный университет Научно-исследовательский институт физики 198904, Россия, Санкт-Петеpбуpг, Петpодвоpец, ул. Ульяновская, 1,к.1 E-mail: com974@mail.ru, tsiganov@mph.phys.spbu.ru

Получено 21 сентября 2005 г.

Рассматривается возможность алгоритмического построения переменных разделения в уравнении Гамильтона—Якоби для обширного класса так называемых L-систем, или систем Бененти, на римановых многообразиях постоянной кривизны. Предложена программная реализация алгоритма построения переменных разделения в системе символьных вычислений Maple, рассмотрен ряд примеров использования данной программы.

Ключевые слова: интегрируемые системы, уравнение Гамильтона—Якоби, разделение переменных.

Yu. A. Grigoryev, A. V. Tsiganov Computing of the separated variables for the Hamilton—Jacobi equation on a computer

We discuss an algorithm for calculating of the separated variables for the Hamilton—Jacobi equation for the wide class of the so-called L-systems on the Riemann manifolds of the constant curvature. We suggest a software implementation of this algorithm in the system of symbolic computations Maple and consider several examples.

Keywords: integrable systems, Hamilton—Jacobi equation, separation of variables.

Mathematical Subject Classifications: 37J35, 37K10, 70H20

1. Введение

Рассмотрим риманово многообразие Q с локальными координатами q = (q1, q2,..., qn) и положительно определенным метрическим тензором G.

На кокасательном расслоении T*Q многообразия Q стандартным образом определим канонические координаты (p, q) и рассмотрим динамическую систему с функцией Гамильтона натурального вида

П

H = T(p,q) + V(q) = J2 (q)PiPj + V(q). (1.1)

i,j=1

Здесь gij (q) — компоненты метрического тензора G, а V(q) — потенциальная энергия системы, которая является гладкой функцией на римановом многообразии Q, каноническим образом поднятой до функции на всем фазовом пространстве T* Q.

Почти в каждом из учебников по классической механике можно прочитать, что одним из самых универсальных методов интегрирования уравнений движения для конечномерных интегрируемых систем классической механики был и остается метод Гамильтона—Якоби — метод разделения переменных. Успех в применении этого метода всегда связан с удачным выбором системы координат, в которых происходит разделение переменных в уравнении Гамильтона—Якоби

H (p,q)= E. (1.2)

Согласно определению Якоби [9], в уравнении Гамильтона—Якоби имеет место разделение переменных,

если существует полный интеграл вида

n

S (Q,ai ,...,an) = ^2Si (Qi ,ai,... ,an), (1.3)

i= 1

где i-тое слагаемое Si зависит только от i-той координаты Qi и от n параметров а1,..., ап. Здесь предполагается, что переменные разделения Q = (Q1,..., Qn) являются координатами Дарбу, т.е. {Qi; Qj} = 0.

В этом случае переменные Q называются переменными разделения, а сопряженные им импульсы находятся из уравнений Якоби

о _ dSi(Qi,aап)

1 dQi ' { ’

Эти уравнения и их более симметричную форму

Фi(Qi,Pi, а1 ,...,an) = 0 (1.5)

называют разделенными уравнениями — дословно “aequatio separata” [7]. Всюду далее под разделением переменных мы будем понимать разделение переменных в уравнении Гамильтона—Якоби, которое позволяет найти полный интеграл этого уравнения вида (1.3).

Главным недостатком метода разделения переменных считается тот факт, что нахождение переменных разделения для каждой конкретной интегрируемой системы является своего рода искусством. Единственным общим результатом, применимым ко всем интегрируемым системам, является следующий критерий, предложенный в работе Леви-Чивита [11]: уравнение Гамильтона—Якоби (1.2) интегрируемо методом разделения переменных, если функция Гамильтона H(P, Q) удовлетворяет следующим n(n — 1)/2 уравнениям при i = j :

dH dH d2H dH dH d2H dH dH d2H dH dH d2H

+ = (1-6)

дР, дР дд,дд, дР, дд, дд,дР дд, дР дР,дд, дд, дд, дР,дР

До недавнего времени исследование системы уравнений Леви-Чивита (1.6) считалось практически невыполнимой задачей в том числе и потому, что сформулирован этот принцип с использованием неизвестных переменных разделения. Однако для некоторых частных классов интегрируемых систем уравнения Леви-Чивита недавно были переписаны в инвариантной геометрической форме, которая не зависит от выбора координатной системы.

Подчеркнем, что единого для всех интегрируемых систем алгоритма построения переменных разделения по-прежнему не существует. Тем не менее, создано несколько алгоритмов решения этой проблемы, претендующих на универсальность и применимых к некоторым достаточно обширным семействам интегрируемых систем [6, 12, 2, 13]. Для этих систем мы уже можем перейти от искусства разделения переменных к практической технологии разделения переменных.

2. Построение переменных разделения с помощью координатных преобразований

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

д = /(<?), р = /'(<?) р, / = с/ъ...,/п). (2.1)

Функции Гамильтона натурального вида (1.1) ковариантны относительно точечных преобразований, т.е. точечные преобразования изменяют функции gгj (д) и V(д), но не меняют форму гамильтониана:

П

н = т(р, д) + V(д) = ]Т ^(д)р.р + V(д). (2.2)

*,,=1

Подставляя функцию Гамильтона (2.2) в уравнения Леви-Чивита (1.6), получим

1 п 1

- р^р.р^-ёгадгё*даёк1 + дгдаё*ёкгё1а - дгёиЕ?гдаёы - д8ё{г^дгёк^ +

^ /1 .. 1 .. .... . , ч (2.3)

+ ^2Р^{2*Гадг*'3даУ + 2ёГЗдгУд^13 + д^уёгг^а ~ дгёгзёзгд8У - К '

+grsдr VдsV = 0.

Здесь дг = д/ддг, д8 = д/дд8 и г = в.

Приравнивая нулю коэффициенты при импульсах Р. в уравнениях (2.3), получим уравнение для метрики gгj и два уравнения, в которые входят gгj и потенциал V. Кроме того, в эти уравнения неявно входит искомая нами функция /(д) (2.1), которая определяет связь исходных физических координат д и переменных разделения д.

В работах Штеккеля [14], Эйзенхарта [6], Калнинса и Миллера [12], Бененти [2] именно этот вариант уравнений Леви-Чивита (1.6) был переписан в инвариантной форме с использованием тензорного анализа на римановых многообразиях.

ЗАМЕЧАНИЕ. В работах [8, 16] установлена тесная взаимосвязь между теорией разделения переменных в уравнении Гамильтона-Якоби в ортогональной системе криволинейных координат на римановых многообразиях и теорией дифференциальных уравнений в частных производных гидродинамического типа. Используя это наблюдение, можно доказать, что нахождение переменных разделения эквивалентно нахождению инвариантов Римана для уравнений гидродинамического типа.

2.1. Тензора Киллинга

При рассмотрении тензоров мы будем применять стандартные соглашения о суммировании и ранге тензоров. Например, производится суммирование от 1 до п по каждому индексу встречающемуся дважды, один раз наверху, другой раз внизу. Кроме этого, будем обозначать одной и той же буквой К = (Ку) = = (К) = (Ку) тензора типа (2,0),(1,1),(0,2), используя для поднятия и опускания индексов метрический тензор.

Рассмотрим риманово пространство постоянной кривизны Q, = п с положительно определен-

ным контравариантным метрическим тензором

■ • д

G = gl:,дi®дj, дк = -^. (2.4)

Любой контравариантный симметрический тензор К на многообразии Q можно взаимнооднозначно отождествить с однородным полиномом в кокасательном расслоении Т*Q по правилу

К = (К*-') ^ Рк = К*"' р ■■■ ру .

Для тензора нулевого ранга, т.е. функции / на многообразии Q, положим Ру = /, где / - каноническое поднятие функции с Q на Т* Q.

Тензором Киллинга К ранга і называется симметричный (і, 0) тензор в пространстве Q, удовлетворяющий тензорному уравнению Киллинга

[К, С]=0 ^ {Рк ,РС} = 0. (2.5)

Согласно определению, Рк является интегралом движения для геодезического потока на римановом многообразии Q.

Конформным тензором Киллинга Ь ранга і называется симметричный (і, 0) тензор в пространстве Q, удовлетворяющий тензорному уравнению

{Рь,Рс} = сРс . (2.6)

Если К — тензор Киллинга, то тензор Ь = К + /(д)С является конформным тензором Киллинга. В этом случае

с = V/ = g' р' д/,

и, поэтому, тензор Ь называют тензором градиентного типа, а / называют потенциалом.

Конформный тензор Ь с взаимно простыми собственными значениями, для которого кручение Ний-енхейса равно нулю

Тт = 2Ь?да' - 2Ь£д<= 0

называют Ь-тензором, или тензором Бененти. Потенциал Ь-тензора равен его следу, т. е. / =<гі = ігасе(Ь).

2.2. Разделение переменных для L-систем

Основные результаты работ [14, 6, 12, 2], можно сформулировать в виде следующей теоремы:

Теорема 1. Уравнение Гамильтона-Якоби (1.2) допускает разделение переменных в ортогональной системе криволинейных координат, если и только если существует тензор Киллинга К второго ранга с простыми собственными значениями и нормальными собственными векторами такой, что

^(К^У) = 0. (2.7)

Здесь d — внешняя производная.

Подчеркнем, что если условие существования тензора Киллинга К с заданными свойствами и уравнение d(KdV) = 0 переписать в координатном виде с использованием переменных разделения, то мы получим исходные уравнения Леви-Чивита (1.6)—(2.3).

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

Существование тензора Киллинга К с заданными свойствами полностью зависит от выбора кинетической энергии системы Т, т.е. от метрики риманового многообразия. Существование одного такого

тензора Киллинга эквивалентно существованию п тензоров Киллинга Кт с простыми собственными значениями и нормальными собственными векторами, которые коммутируют как линейные операторы и находятся в инволюции относительно соответствующих скобок Пуассона. Линейное пространство, образованное тензорами Кт, называют пространством Киллинга-Штеккеля.

В 1992 году Бененти предложил рекуррентную процедуру построения одного частного семейства пространств Киллинга-Штеккеля [2].

Теорема 2. Если Ь является L-тензором на римановом многообразии Q, то п тензоров

являются (2,0) тензорами Киллинга второго ранга с простыми собственными значениями и нормальными собственными векторами, которые попарно коммутируют друг с другом.

Входящие в определение (2.8) функции стто являются симметрическими полиномами степени т от собственных значений тензора Ь, то есть они равны коэффициентам характеристического полинома

Уравнение ^(К^У) = 0 обычно рассматривается как уравнение для нахождения потенциалов V, которые можно добавить к кинетической энергии Т на данном римановом многообразии с сохранением свойства разделимости переменных [6, 12, 13]. Легко доказать, что из уравнения ^(К^У) = 0 следует, что ^(КТО^У) = 0 для любого т. Поэтому, согласно [2], можно ввести семейство потенциалов Ут, которые являются решениями уравнений

которые образуют семейство интегралов движения в инволюции для заданного гамильтониана Н (1.1).

Согласно [10], интегралы движения Нт (2.10) также можно найти в виде решений рекуррентных уравнений

Здесь N - оператор рекурсии, который является каноническим поднятием тензора Ь на кокасательное расслоение Т* Q по правилу

Построение тензоров Киллинга Кт и соответствующих интегралов движения Нт тесно связано с разделением переменных согласно следующей теореме.

Теорема 3. Собственные значения ^ тензора Ь являются переменными разделения для интегрируемой системы с интегралами движения Нт.

Интегрируемые системы, допускающие разделение в этих переменных, называют Ь-системами, системами Бененти, или системами кофакторного типа.

ЗАМЕЧАНИЕ. Подчеркнем, что L-системы не исчерпывают всего множества систем Штеккеля, допускающих разделение переменных в ортогональной системе криволинейных координат на римановом многообразии 2. Например, конические координаты не являются нормальными координатами для какого-то либо тензора Ь. Для рассмотрения подобных координатных систем используют понятие пучка L-тензоров [4, 13].

т

(2.8)

П

ае^АІ — Ь) — ^2 ^тА”-т

т=0

(2.9)

и определить интегралы движения

П

(2.10)

^Нт+1 — ^¿Нт + От+1 ¿Н, т — 1,...,п — 1, Н„ = 0 .

(2.11)

(2.12)

ЗАМЕЧАНИЕ. Согласно [15], (1,1)-тензор — цО) является полиномом степени п — 1 по коэффи-

циенты которого с точностью до знака совпадают с тензорами Кт (2.8). Напомним, что то^А) определен следующим образом

то^А)А = Acof(A) = ііе1;(А)С.

В силу этого системы с пространством Киллинга-Штеккеля вида (2.8) и называют системами кофакторного типа.

2.3. Алгоритм построения переменных разделения

Итак, для заданного гамильтониана Н (1.1) в рассмотренном выше методе построения переменных разделения необходимо найти тензор Киллинга Ь и нетривиальные решения уравнений 1(К1У) — 0 или (2.9)—(2.11).

Для создания алгоритма построения переменных разделения условие существования тензора Кил-линга Ь необходимо также представить в виде условия существования нетривиального решения некоторого уравнения. Искомое уравнение было построено в работе [5], в которой доказана следующая теорема.

Теорема 4. Тензор Ь является симметрическим тензором Киллинга градиентного типа c нулевым кручением и простыми собственными значениями относительно метрического тензора С, если и только если

1(£хт в — Т!^) — 0. (2.13)

Здесь £ — производная Ли вдоль геодезического векторного поля Хт, ст1 — 1хЬ — симметрический полином первой степени и

П

в — Е Ь Рі^ (2.14)

¿,¿=1

является L-деформацией стандартной 1-формы Лиувилля в0 — ^ р-1д', не зависящей от выбора координатной системы.

Подставив тензор Киллинга К1 (2.8) в уравнение 1(К1У) — 0, мы получим уравнение, которое по своей структуре аналогично уравнению (2.13)

1(£х^ в — ) —0. (2.15)

Доказательство может быть найдено, например, в работе [1].

Для того, чтобы уменьшить количество промежуточных вычислений, мы воспользуемся следующим выражением для производной Ли £ вдоль векторного поля X

£х — іх 1 + 1 іх .

Здесь іх — оператор крюка и 1 — внешняя производная. Подставим это выражение в уравнения (2.132.15). Так как 12 — 0, то

1£х — 1 іх 1 + 12 іх — 1 іх 1, и уравнения (2.13) и (2.15) имеют вид

1(іхт 1в — Т1ст1) — 0, (2.16)

1(іх^ 1в — У1ст1) — 0. (2.17)

Итак, алгоритм построения переменных разделения Q для данной нам Ь-системы сводится к нахождению нетривиальных решений уравнений (2.16) и (2.17) относительно функций Ь'(д) на римановом многообразии Q и нахождению собственных значений тензора Ь. В этом алгоритме для построения переменных разделения, вместо неизвестного вектора f (2.1) на римановом многообразии Q, который неявно входит в исходные уравнения (2.3), мы будем искать тензор Киллинга Ь второго рода со специальными свойствами.

В следующем разделе мы представим реализацию данного алгоритма для построения Ь-тензоров и соответствующих переменных разделения для Ь-систем.

ЗАМЕЧАНИЕ. Если для данного гамильтониана уравнения (2.16-2.17) не имеют решения или решение тривиально, то это означает, что данная интегрируемая система не является L-системой. Отсутствие решения вовсе не означает, что система не допускает разделения переменных.

ЗАМЕЧАНИЕ. Вторая дифференциальная 1-форма в естественным образом порождает вторую пуассонову структуру на многообразии T*Q [10, 1]. Соответствующий оператор рекурсии N (2.12) является каноническим поднятием тензора L на кокасательное расслоение T* Q. В силу этого предложенный алгоритм можно рассматривать как часть более общего алгоритма нахождения переменных разделения для би-гамильтоновых интегрируемых систем.

3. Программа для построения переменных разделения

Далее приведена программная реализация описанного выше алгоритма построения переменных разделения, выполненная в среде символьных вычислений Maple v.9.5. Строки программы помечены символом >, остальное представляет собой развернутый комментарий. В виде исполняемого Maple-файла данная программа размещена для свободного доступа по адресу

http://www.maplesoft.com/applications/app_center_view.aspx?AID=1686

В данной программной реализации нами использовался пакет liesymm, который является стандартной встроенной библиотекой системы символьных вычислений Maple. Он ориентирован на построение и работу с дифференциальными формами, получающимися из известных дифференциальных уравнений в частных производных. В данном разделе нами осуществляется обратная процедура: мы построим дифференциальные уравнения, отвечающие условиям (2.16-2.17), записанным в терминах дифференциальных форм.

Итак, загрузим необходимую нам библиотеку встроенных процедур и зададим размерность конфигурационного пространства Q:

> with(liesymm):

> n:= 2;

n := 2

Координаты на римановом многообразии Q будем обозначать q = (q1,..., qn), а сопряженные им импульсы p = (pi ,...,p„):

> q:=seq(q|І і,i=l..n): p:=seq(p||і,i=l..n): var:=q,p:

> setup(var);

[q1, q2, p1, p2]

В уравнения (2.13-2.15) входят векторные поля, части функции Гамильтона и компоненты конформного тензора Киллинга L. Начнем с определения следа ст1 тензора L и соответствующей 1-формы в (2.14)

> sigma:=add(L[і,і](q),і=1..п);

а := Li, i(q1, q2) + ¿2,2 (q1, q2)

> theta:=add(add(L[i,j](q)*p|Ii*d(q|Ij),i=l..n),j=l..n);

в := Li, i(q1, q2) p1 d(q1) + ¿2, i(q1, q2) p2d(q1)

+Li, 2 (q1, q2) p1 d(q2) + ¿2,2(q1, q2) p2d(q2)

Здесь Lj (q) - искомые нами функции на римановом многообразии Q.

Далее, для нахождения интегралов движения с помощью рекуррентных соотношений (2.11) нам потребуется тензорное поле N типа (1,1), действие которого на произвольную 1-форму определено формулой (2.12):

> Nq:={seq(d(q[k])=add(L[k,j](q)*d(q[j]) ,j=l..n),k=l..n)}:

> Np:={seq(d(p[k])=add(L[j,k](q)*d(p[j])

> -pI Ij*add((diff(L[j,i](q),q[k])-diff(L[j,k](q),q[i]))*d(qI Ii),i=l..n),

> j=l..n),k=l..n)>:

На втором шаге вводятся канонический тензор Пуассона 3?, который необходим нам для построения векторных полей

> ed:=array(identity, l..n,l..n): u:=array(sparse,1..n,l..n):

> P:=linalg[blockmatrix](2,2,[u,ed,-ed,u]);

P :=

0 0 10

0 0 0 1

-1 0 0 0

0 -10 0

и канонические скобки Пуассона

> PB:=proc(f,g) options operator, arrow:

> add( diff(f,p||i)*diff(g,q||i)- diff(f,qIIi)*diff(g,p||i), i=l..n)

> end:

Используя эти определения, мы теперь можем вычислить векторное поле XT = 3?dT(p, q)

> dT:=array(1..2*n):

> for i from 1 to 2*n do dT[i]:=diff(T(var),var[i]): end do:

> X_T:=evalm(P&*dT);

XT := 7-—— T(ql, q2, pi, p2), ^r~ T(ql, q2, pi, p2),

L dpi dp2

~ T(qi, q2, pi, p2)), P1 > P2))

и векторное поле XV = 0PdV(q)

> dV:=array(1..2*n):

> for і from 1 to 2*n do dV[i]:=diff(V(q),var[i]): end do:

> X_V:=evalm(P&*dV);

XV :=

°. o. -(4lv(9l’92))’"(^v(9l’92))

Теперь у нас все готово для того, чтобы определить уравнения (2.16-2.17) и построить с помощью них искомую систему дифференциальных уравнений в частных производных.

Начнем с уравнения (2.17), которое по объему вычислений требует значительно меньше ресурсов. Подставляя векторное поле Ху и 1-форму в в выражение і(Х)

> idT: =wcollect (Ь.оокСсК'ЬЬ.еЪа) ,Х_У) ) :

гсІТ := ( - Ь2,і(д1, ?2) У(?1, q2))(ql, q2, р1, р2) -

~ д2) |?2))^1, Я2, Р!. р2))с%1) +

+ ( - Ь2,2(ді, q‘i) (-т^р у(«1> ч2, рі, р2) -

~ q2) (щ- У^1, q2))(ql, q2, р1, р2))с%2),

мы получим соотношение, в которое входят производные, формально зависящие от всех переменных фазового пространства. В действительности, эти производные зависят только от координат на конфигурационном пространстве. Это особенность реализации встроенной процедуры i(X), которую нам необходимо исправить.

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

> Trans:={seq(diff(V(q),qlIi)(var)=diff(V(q),qlIi),i=l..n)}:

> idT:=subs(Trans,idT);

idT:= ( - L2a(ql,q2)^V(ql,q2)-LM(ql,q2)^IV(ql,q2))d(ql)

+ ( - Ь2,2(gl, q2) V(ql, q2) - Lii2(ql, q2) ?2))d(?2)

Приравнивая нулю коэффициенты при всех независимых 2-формах в (2.17), мы получим следующую систему алгебраических и дифференциальных уравнений:

> SysEq:=anmil(d(idT-V(q)*d(sigma)), [var] ): nops(SysEq);

При n = 2 эта система состоит из одного уравнения, при n = 3 система состоит из трех уравнений и т.д. Полученные уравнения мы поместим в подготовленный нами пустой список уравнений

> ListEq:=NULL:

> for i from 1 to nops(SysEq) do

> ListEq:=ListEq,lhs(SysEq[i]);

> end do :

Теперь перейдём к оставшемуся уравнению (2.16). Как и раньше, приходится вручную упрощать это уравнение после применения встроенной процедуры :

> Eq:=wcollect(d(hook(d(thêta),X_T)-T(var)*d(sigma))):

> Trans:={seq( diff(T(var),var[i])(var)=diff(T(var),var[i]),i=l..2*n)}:

Приравнивая теперь нулю коэффициенты при всех независимых 2-формах в уравнении (2.16), мы получим вторую систему алгебраических и дифференциальных уравнений

> SysEq:=anmil(subs(Trans,Eq),[var]): nops(SysEq);

Уравнение (2.16) при n = 2,3 порождает шесть и пятнадцать уравнений, соответственно. Добавим эти уравнения к ранее подготовленному списку

> for i from 1 to nops(SysEq) do ListEq:=ListEq,lhs(SysEq[i]); end do:

> NumEq:=nops([ListEq]);

NumEq := 7

При n = 2,3 полная система уравнений состоит из семи и восемнадцати уравнений для четырех и девяти функций Lj (q), соответственно. Эти уравнения содержат координаты qj, функции от координат Lj (q), гамильтониан H и импульсы pj, которые входят в определение гамильтониана.

Поэтому на третьем шаге в полученную нами систему уравнений необходимо подставить функцию Гамильтона H для исследуемой нами интегрируемой системы. В качестве первого примера рассмотрим интегрируемую модель Энона—Эйлеса с гамильтонианом

H = р\ + р2 + ^(qi + q2) + + 2ql. (3-1)

Соответствующая кинетическая и потенциальная энергии T и V определяются следующим образом:

> Tn:=add(p||i~2,i=l..n);

> Vn : =1/2* (ql~2+q2'-2)+ql'-2*q2+2*q2~3 ;

Подставляя данные части функции Гамильтона

> Н:={ T(var)=Tn, V(z)=Vn };

в полученную нами ранее систему уравнений ListEq и приравнивая нулю коэффициенты при степенях переменных pj, мы получим искомую систему алгебраических и дифференциальных уравнений в частных производных System для элементов тензора Lj (z1;..., zN _i):

> System:= NULL:

> for i from 1 to NumEq do

> u:=expand(dvalue(subs({ T(var)=Tn, V(z)=Vn },ListEq[i]))):

> System:=System,coeffs(u,{p}):

> end do:

> NumSys:=nops({System});

NumSys := 13

При n = 2,3,4 результирующая система уравнений состоит из 13, 51 и 136 уравнений, но не все эти уравнения линейно независимы.

ЗАМЕЧАНИЕ. Для произвольного потенциала V (q) конечная система уравнений обычно сильно переопределена и имеет только тривиальное решение, что означает, что переменные разделения для данной системы не могут быть построены с помощью точечных или координатных преобразований для L-систем.

В нашем случае система уравнений System легко решается с помощью стандартной процедуры pdsolve

> âns:=pdsolve( {System}, {seq(seq(L[i,j](z), i=l..n),j=l..n)} );

Соответствующее решение

Ans := {Llyl(ql, q2) = + C2, Liy2(ql, q2) = —

L2,2 (ql, q2) = Ci q2 + C2, L2jl(ql, q2) = —y“}

зависит от двух произвольных постоянных Ci,2. Подставляя это решение в L-тензор, мы получим

> Trans :={seq(seq(L[i,j](q)=L[i,j],i=l..n),j=l..n)}:

> ins:=simplify(map2(subs,Trans,ins)):

> L:=array(1..n,1..n):

> L:=map2(subs,ins,evalm(L));

L :=

Cl,‘ ClqZ + Ch

2

Собственные значения этого Ь-тензора по определению являются переменными разделения

> СЗ: =Ипа1§ [е1§еп¥а1иез] (Ln) ;

„ 3Сх _ С\ Ч2 \/9 С\ - 24 С\ q2 + 16 С\ q22 + 16 С\ ql2

Ч '■= —Б-----ь с2 н--------ь

8 z 2 8 3Ci Ciq2 фС1- 24 Cl q2 + 16 Cl q2 2 + 16 Cl ql 2

+ c2 + -

8

2

8

Эти координаты совпадают со сдвинутыми параболическими координатами, которые действительно являются переменными разделения для системы Энона-Эйлеса [13]. Построение этих переменных с помощью компьютера заняло менее одной минуты машинного времени.

Строго говоря, мы построили целое семейство эквивалентных систем переменных разделения, отличающихся произвольными постоянными Cij2. Следует напомнить, что две системы переменных разделения называются эквивалентными, если соответствующие решения уравнения Гамильтона-Якоби порождают одно и то же лагранжево расслоение фазового пространства T*Q.

Теперь обсудим построение соответствующих интегралов движения с помощью семейства конформных тензоров Киллинга (2.8). Определим полиномы стто

> for i from 0 to n do

> sigma[i]:=coeff(linalgEdet](lambda*ed-L).lambda,n-i):

> end do:

и базис из тензоров Киллинга Km (2.8)

> for m from 1 to n-1 do

> K||m:=evalm(add( sigma[m-k]*L~k,k=0..m)):

> end do:

Напомним, что в нашем случае gij = 1 и, таким образом, Lij = Lj.

Кроме тензоров Km уравнение (2.9) содержит внешние производные от неизвестных нам пока потенциалов , которые также необходимо определить

> for m from О to n-1 do

> dV||m:=array(l..n): for i from 1 to n do

> dVl|m[i]:=diff(VI Im(q),q[i]):

> end do:

> end do:

> dV0:=map2(subs,V0(q)=Vn,dV0):

Теперь у нас все готово для того, чтобы построить и решить уравнения (2.9), что и позволит нам определить дополнительные интегралы движения (2.10)

>for m from 1 to n-1 do

> eqV:=evalm(dVNm-K||m&*dVO):

> Ans:=pdsolve({seq(eqV[i], i=l..n)},V||m(q)):

> H[m]:=add(add( K||m[i,j]*p||i*p|Ij ,i=l..n),j=l..n)+subs(Ans,V||m(q)):

> end do:

Мы не будем выписывать полученные интегралы движения Н [т] явно, но проверим, что эти функции действительно находятся в инволюции с гамильтонианом

> Н[0]:=Tn+Vn;

Н0 := pi2 + Р22 + ^ ^ + ql2 q2 + 2 q23

Действительно, все скобки Пуассона между Н СО] и Н [т] равны нулю

> for т from 1 to n-1 do

> ZERO:=simplifу(PB(H[0],H[m]));

> end do;

ZERO := 0

ЗАМЕЧАНИЕ. Если число степеней свободы не очень велико, например n ^ 10, то для полиномиальных потенциалов переопределенную систему уравнений System можно решить и с помощью стандартной процедуры pdsolve на обычном персональном компьютере за разумное время. При n > 10 необходимо использовать специальные компьютеры или специальные системы символьных вычислений для решения полученных нами сильно переопределенных систем уравнений.

4. Примеры

4.1. Ангармонический осциллятор

Рассмотрим n-мерное евклидово пространство Q = Rn с декартовыми координатами qг такое, что метрика равна ds2 = gjdq*dqj = Y1 dq*2. Гамильтониан ангармонического осциллятора имеет вид

n n / n \ 2

H = Ep2 + E a* q*2 + E q*2 . (4.1)

¿=i *=1 V *=1 )

В системе Maple соответствующая кинетическая и потенциальная энергии определяются следующим образом:

> Tn:=add(pИi~2,i=l..n);

> Vn:=add(a||i*q||i~2,i=l..n)+(add(q||i~2,i=l..n))~2;

Как и для системы Энона-Эйлеса, нам необходимо подставить этот гамильтониан в систему уравнений и затем решить переопределенную систему уравнений , в которую уже не входят момен-

ты pj. Соответствующее решение имеет вид

Lj (q) = (Ci + a* — an)^jj + C2 q*qj

При C1 = an и C2 = 1 собственные значения Q * соответствующего тензора L определяются уравнением

det(L - Л) _ vp qj2 =ttA-Q8

n?=i(A-ai)“

которое совпадает с известным определением эллиптических координат в евклидовом пространстве Rn.

Перейдем теперь к построению интегралов движения. В отличие от системы Энона-Эйлеса, для построения интегралов мы воспользуемся методами би-гамильтоновой геометрии. Начнем с того, что подставим найденный нами тензор L в общее определение (2.12)

> Nqn:=wcollect(simplify(map2(subs,Ans,Nq))):

> Npn:=wcollect(simplify(map2(subs,Ans,Np))):

Затем вычислим симметричные полиномы стто от собственных значений тензора L:

> for m from 0 to n do

> sigma[m]:=coeff(linalgEdet](Ln-lambda*ed),lambda,n-m):

> end do:

Теперь у нас все готово для того, чтобы построить и решить рекуррентные соотношения (2.11)

> unassign(H): H[n+1]:=0: H[0]:=Tn+Vn;

> for m from n-1 by -1 to 1 do

> Eq:=wcollect(dvalue( d(H[m+l])-

> subs(Nqn union Npn,d(H[m](var)))-sigma[m+1]*d(H[0]))):

> Sys:=annul(Eq,[var]):

> Ansi:=pdsolve(Sys,H[m](var));

> H[m]:=subs(Ansi,H[m](var)):

> end do:

Мы не будем выписывать полученные интегралы движения Н [т] явно, но проверим, что эти функции действительно находятся в инволюции с гамильтонианом H (4.1)

> for m from 1 to n-1 do

> ZERO:=simplify(PB(H[0],H[m]));

> end do:

ZERO := 0

Все скобки Пуассона между Н СО] и Н [т] равны нулю.

4.2. Эйлеровская задача двух центров

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

H = pi + р2 -

й!

+ ■

«2

y?{qi ~c)2 + ql y/(qi +с)2 +qj

В отличие от предыдущих примеров в этом случае потенциал V(д) является алгебраической, а не полиномиальной функцией от координат.

Решение переопределенной системы уравнений Зу81ет имеет вид

L

(q2 - c2)C2 + Ci 5152^2

qiq2 C2 qf C2 + Ci

Если С1 = 1/2 с2 и С2 = 1, то собственные значения матрицы Ь являются стандартными эллиптическими

координатами

1 -

q?

«2

_ (^ — — Я2)

Л + с2/2 Л — с2/2 " (Л — с2/2)(Л + с2/2) '

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

4.3. Цепочка Тоды

Пусть конфигурационное пространство Q совпадает с п-мерным евклидовым пространством К”. В декартовых координатах функция Гамильтона для периодической решетки Тоды имеет вид

Н — \ + y~!exp(g» — Чг+і), q-n+i = qi ■

(4.3)

i= 1 i= 1

При n = 2 существует единственное решение полной системы уравнений System

L

Ci C2 C2 Ci

Конечно, собственные значения этой матрицы не являются переменными разделения. Тем не менее, в этом случае число свободных параметров равно размерности риманова многообразия Q и этот факт позволяет нам все равно построить переменные разделения. Для этого необходимо использовать идею близкого к нашему алгоритма построения переменных разделения, подробно описанного в работе [13]. Действительно, приведем матрицу Ь к диагональному виду с помощью преобразования подобия

L = V-idiag(Ci + C2,Ci - Ci)V, V:

11 -1 1

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

Q = Vq, ^ Qi = qi + q2, При n = 3 решение системы уравнений System имеет вид

Q2 = qi + q2.

L-

C1 C2 C2 C2 C2 C22 C2 C2 C1

Число свободных параметров меньше размерности риманова многообразия и, поэтому, такое решение называется тривиальным [13]. Если решение тривиально, это означает, что для данного гамильтониана невозможно построить переменные разделения с помощью преобразований конфигурационного пространства Q.

Тем не менее, необходимо подчеркнуть, что для цепочки Тоды уравнение Гамильтона-Якоби допускает разделение переменных, но для построения переменных разделения необходимо использовать более широкий класс канонических преобразований фазового пространства T* Q.

4.4. Система Неймана

Конфигурационное пространство Q для системы Неймана является сферой Sn единичного радиуса в евклидовом пространстве R"+1. С помощью декартовых координат х = (xi,..., x„_|_i) в R”+1, риманово многообразие §" определяется условием \х\ = у/(ж, х) = 1.

Система Неймана описывает частицу, движущуюся по поверхности сферы Sn в поле квадратичного

потенциала V = ^ J2 aixi- Нам будет удобно производить все вычисления в следующих локальных координатах на сфере: q» = ж2, i = 1,..., n. В этих координатах гамильтониан системы Неймана H = T + V имеет вид

n n n

1

2

i=1 i<j i=1

Здесь p — импульсы, канонически сопряженные координатам qj, а a1,..., «n+1 — произвольные числовые параметры, входящие в определение потенциала. В системе Maple этот гамильтониан можно задать с помощью следующих команд

> Tn:=2*add(q||i*(l-q||i)*p||i~2,i=l..n)

> -4*add(add(q| |i*q|Ij*pI Ii*pI Ij,j=1+1..n),i=l..n);

> Vn:=l/2*add((a||i-a|I(N+l))*q||i,i=l..n);

Если n = 2, то решение переопределенной системы уравнений System равно

„ 91 («3 — «1) + «1 — «2 „ „

С2----------------------- - -Ci qi Ci

«2 - «3

q2 («1 - «з) ^ ^ ^

—----------Cl q2 Cl + C2

«2 - «3

T = 2 ^2 qi(l - qi)pi ~ 4 ^ qiQjPiPj, V = - ^(ai ~ ап+і)Чі ■

L

Напомним, что Ь* (д) — компоненты тензора Ь, один раз ковариантного и один раз контравариантно-го. Для того, чтобы поднять или опустить индекс, необходимо использовать метрический тензор С. Для неплоских метрик это приводит к тому, что Ь-матрицы, полученные в результате работы программы, несимметричны, хотя и отвечают симметричному тензору Киллинга.

Собственные значения ^12 полученного Ь-тензора являются нулями характеристического полинома, который при С1 = а3 — а2 и С2 = а2 имеет вид

— Л) ді д2 1 — ді — д2

ел = —--------------= ---------Ь -----------г-------.

Пд=і (А — аі) Л — аі А — а2 А — аз

В исходных избыточных координатах х* в К”+1 переменные разделения определяются следующим образом:

^ _1Е., <*-<*> *

Л - а,- ГГ+1 (Л _ а л’ ^

3= Л _ а, З (Л _ а,)

з

1^=1 ™з)

Это известное определение сфероконических координат или эллиптических координат на сфере §”.

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

> Nqn:=wcollect(simplify(map2(subs,Ans,Nq))):

> Npn:=wcollect(simplify(map2(subs,Ans,Np))):

> for m from 0 to n do

> sigmaEm]:=coeff(linalgEdet](Ln-lambda*ed),lambda,n-m):

> end do:

> unassign(H): HEn+l]:=0: HEO]:=Tn+Vn;

> for m from n-1 by -1 to 1 do

> Eq:=wcollect(dvalue( d(HEm+l])-

> subs(Nqn union Npn,d(HEm](var)))-sigmaEm+l]*d(HE0]))):

> Sys:=annul(Eq,Evar]):

> Ansi:=pdsolve(Sys,HEm](var));

> HEm]:=subs(Ansl,HEm](var)):

> end do:

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

> for m from 1 to n-1 do

> ZERO:=simplify(PB(HEO],HEm]));

> end do;

ZERO := 0

В данном примере, когда риманово многообразие снабжено неплоской метрикой, все вычисления на обычном настольном компьютере при n = 2,3 занимают примерно минуту и восемь минут, соответственно.

Замечание: В отличие от вычислений, приведенных в работе [1], мы непосредственно решаем уравнения (2.13-2.15), а не применяем для решения подстановку специального аффинного вида. Это позволяет говорить о единственности полученного решения.

4.5. Модель Якоби-Калоджеро

Пусть конфигурационное пространство Q — трехмерное евклидово пространство R3 с декартовыми координатами q = (q1; q2, q3). Функция Гамильтона равна

3

2

н = ^3p2 + (q1 - q2) 2 + (q2 - 9з) 2 + (9з - 91)"

j=1

В случае п = 3 данный гамильтониан рассматривался Якоби, если п > 3, то данную систему принято называть системой Калоджеро.

Решение полной системы уравнений ЭузЪеш зависит от четырех произвольных параметров

L

92C1 + 2q1 C2 + C3 9192C1 + (91 + 92)C2 + C4 919361 + (91 + 93)62 + C4

9192 C1 + (91 + 92)62 + C4 q| C1 +29262 + C3 92 9361 + (92 + 93)62 + C4

91 9з61 + (91 + 93)62 + 64 929з61 + (92 + 93)^2 + 64 q261 + 29362 + 63

Собственные значения этого Ь-тензора являются переменными разделения в соответствующем уравнении Гамильтона—Якоби

> Q:=linalgEeigenvalues](Ь);

Число свободных параметров С больше, чем размерность риманова многообразия. Это означает, что данная интегрируемая система является вырожденной или суперинтегрируемой системой, которая допускает разделение переменных в нескольких различных координатных системах. Действительно, при различном выборе свободных параметров Сі,..., С4 можно получить собственные значения тензора Ь, которые совпадут с сплющенными сфероидальными, вытянутыми сфероидальными, сферическими и пара-болоидальные координатами, также можно получить координаты параболического цилиндра. Детальное описание этих координатных систем и их применение к интегрированию уравнений движения в модели Якоби—Калоджеро может быть найдено в работах [13, 3].

4.6. Рациональный потенциал

Рассмотрим движение точки в обычном евклидовом пространстве К” с декартовыми координатами 0 = (оі, 02, • • •, 0”). Функция Гамильтона равна

В отличие от предыдущих примеров потенциал V(о) является рациональной, а не полиномиальной функцией координат.

Сингулярная для данного потенциала поверхность является эллипсоидом, который задается уравнением

При с > 0 все траектории, начинающиеся внутри данного эллипсоида, остаются внутри навсегда.

На обычном настольном компьютере при n > 2 стандартная процедура pdsolve не может решить систему уравнений System за разумное время, т.е. решение занимает более нескольких десятков минут. Это увеличение времени работы связано с общими особенностями системы Maple, которые всегда возникают при работе с корнями, рациональными степенями и комплексными функциями.

Эти системные особенности легко обойти, если сначала найти общее решение уравнения (2.16), которое описывает геодезическое решение и не содержит рациональных функций, входящих в определение потенциала. Конечно, затем это общее решение необходимо подставить в уравнение (2.17) и найти решение, описывающее тензор L для данного конкретного потенциала. В рассматриваемом нами примере такое изменение алгоритма решения позволяет сократить затраты машинного времени до нескольких минут при n ^ 10.

Для евклидова пространства общее решение системы уравнений, полученных из уравнения (2.16) и описывающих движение по геодезическим, имеет вид

Здесь коэффициенты а, Ау, Ду и Су — произвольные постоянные, такие-,- что Ау = Ду и Су = Су.

Подставляя это общее решение в систему уравнений, полученную из второго уравнения (2.17), мы получим систему алгебраических уравнений на коэффициенты Ау, Ду и Су, вид которой полностью определяется конкретным потенциалом.

Для рационального потенциала (4.4) стандартная Мар1е процедура зоГуе с легкостью решает эту систему алгебраических уравнений за несколько секунд. При этом собственные значения соответствующего тензора Ь являются стандартными эллиптическими координатами (4.2) в евклидовом пространстве.

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

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

П

(4.4)

n ^2

Программная реализация данного алгоритма осуществлена в системе символьных вычислений Maple. Программа и примеры ее использования, перечисленные в данной работе, размещены-а- для свободного доступа по адресу

http://www.maplesoft.com/applications/app_center_view.aspx?AID=1686

В дальнейшем мы собираемся внести ряд изменений в данную программу и создать библиотеку программ для работы с L-системами. В частности, добавить программы, необходимые для построения и решения соответствующих разделенных уравнений, программы для построения оператора рекурсии, матриц Лакса, преобразований Бэклунда и функций Бейкера-Ахиезера.

Список литературы

[1] Bartocci C., Falqui G., Pedroni M. A geometric approach to the separability of the Neumann-Rosochatius system // Differential Geom. Appl., 2004, V. 21, № 3, p. 349-360.

[2] Benenti S. Orthogonal separable dynamical systems. In: Proceedings of the 5th International Conference on Differential Geometry and Its Applications, Silesian University at Opava, August 24-28, 1992 (O. Kowalski and D. Krupka Eds.)//Differential Geometry and Its Applications, 1993, V 1,p. 163-184.

[3] Benenti S., Chanu C., Rastelli G. The super-separability of the three-body inverse-square Calogero system // J. Math. Phys., 2000, V. 41, p. 465-4678.

[4] Benenti S., Chanu C., Rastelli G. Remarks on the connection between the additive separation of the Hamilton-Jacobi equation and the multiplicative separation of the Schrodinger equation // J. Math. Phys., 2002, V. 43, p. 5183-5253.

[5] Crampin M., Sarlet W, Thompson G. Bi-differential calculi, bi-Hamiltonian systems and conformal Killing tensors // J. Phys. A, 2000, V. 33, p. 8755-8770.

Crampin M. Conformal Killing tensors with vanishing torsion and the separation of variables in the

Hamilton-Jacobi equation // Diff. Geom. Appl., 2003, V. 18, p. 87-102.

[6] Eisenhart L. P. Separable systems of St&ckel // Ann. Math., 1934, V. 35, p. 284-305.

[7] Euler L. Calculi integralis // Ac. Sc. Petropoli, 1768, V. 1-3. (Russian translation: GITL, Moskow, 1956.)

[8] FerapontovE. V., Fordy A. P. Separable Hamiltonians and integrable systems of hydrodynamic type // J. Geom. Phys., 1997, V 21, p. 169-182.

[9] Jacobi C. G. J. Vorlesungen uber Dynamik // Georg Reimer, Berlin, 1866. (Jacobi's lectures on dynamics given in Konigsberg 1842-1843, published by A. Clebsch.)

[10] Ibort A., Magri F., Marmo G. Bihamiltonian structures and Stackel separability // J. Geometry and Physics, 2000, V 33, p. 210-228.

[11] Levi-Civita T. Integrazione delle equazione di Hamilton-Jacobi per separazione di variabili // Math. Ann., 1904, V. 24, p. 383/2397.

[12] Kalnins E.G., Miller W., Jr. Killing tensors and variable separation for Hamilton-Jacobi and Helmholtz equations // SIAM J. Math. Anal., 1980, V. 11, p. 1011-1026.

[13] Rauch-Wojciechowski S, Waksjo C. How to find separation coordinates for the Hamilton-Jacobi equation: a criterion of separability for natural Hamiltonian systems // Math. Phys. Anal. Geom., 2003, V. 6, №4, p. 301-348.

[14] Stackel P. 0ber die Integration der Hamilton—Jacobischen Differential Gleichung mittelst Separation der Variabeln // Habilitationsschrift, Halle, 1891.

[15] Schouten J. A. Ricci Calculus // Springer, Berlin, 1954.

[16] Zakharov V. E. Description of the n-orthogonal curvilinear coordinate systems and hamiltonian integrable systems of the hydrodynamic type. Part I. Integration of the Lame' equation // Duke Math. J., 1998, V. 94, p. 103-139.