Реализация мультискобки Кругликова-Лычагина в системе компьютерной алгебры Maple Текст научной статьи по специальности «Математика»

УДК 517.951 ББК 22.161

РЕАЛИЗАЦИЯ МУЛЬТИСКОБКИ КРУГЛИКОВА-ЛЫЧАГИНА В СИСТЕМЕ КОМПЬЮТЕРНОЙ АЛГЕБРЫ MAPLE

Тычков С. Н.1

(Учреждение Российской академии наук Институт проблем управления им. В.А. Трапезникова РАН, Москва)

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

Ключевые слова: дифференциальные уравнения, мультискобка Кругликова-Лычагина, Maple.

Введение

В работе предлагается реализация мультискобки Кругликова-Лычагина, которая используется для решения вопроса о формальной интегрируемости системы дифференциальных уравнений [1]. Такие задачи часто возникают в дифференциальной геометрии, теории дифференциальных уравнений и геометрической теории управления.

Рассмотрим систему m + 1 дифференциального уравнения в частных производных порядка s

где х = (х1, х2,..., хп) — вектор размерности п; и = (и1, и2,..., ит) — вектор-функция от х; а = (а1,..., ап) — муль-тииндексы длины |а| = а1 + ... + ап. В этих обозначениях под

1 Сергей Николаевич Тычков, аспирант, ([email protected]).

(1)

i = 1,..., m + 1,

д|а|иг/дха мы

понимаем частную производную

д^и

дх^1 дх°2 ...дхПп

Для упрощения записи введем обозначение:

к = 1,..., т.

Ра =

дМ ик

дх0

Под оператором полной производной мы понимаем оператор

ё д ^ к д

/ 1

к,а

дР

где а + 1j — мультииндекс, полученный из мультииндекса а увеличением его ] -й компоненты на единицу.

Рассмотрим дифференциальное уравнение порядка в вида

(2) F (х,и,р\ 1 ,...,Р™ ) =0.

Определение 1. Линеаризацией [2] уравнения (2) называется вектор дифференциальных операторо)

(3)

где

Ьк

т = (ь1

дрО х'

,Ьт) ,

к = 1,..., т,

- = -о ... о-,

х ёа1 х1 хп

о — композиция дифференциальных операторов, a

--, ] = 1,..., п

xj

— оператор полной производной порядка аj по переменной xj•.

Линеаризуя каждое уравнение системы (1), получим линеаризованную систему в виде: /

(4)

ь1

ь1

\

т 1

^т 1

"т+1

Т 2

т

Т 2

ьт+1

т т ьт

и1

ит

0

0

1 2 т

\ ьт+1 ьт+1 . . . ьт+1 /

где Ь — ^-я компонента вектора

Пусть Т>г = (р1,..., , г = 1,..., в, где р — линейный дифференциальный оператор. Множество линейных дифференциальных операторов образует некоммутативную алгебру относительно композиции. 128

Определение 2. Некоммутативным определителем

(Р1,...,Vm) называется любое (выбранное и зафиксированное) обобщение определителя для алгебры дифференциальных операторов:

/ V г>2 ... ^т \

Ф1,..., Vm) =Ndet

I V1 V2 Т>т

\ Vm • • Vm /

,1 т»т

а(1) ... Va(m),

£ (-1)а^(1)...Vm

а£Яп

где Бт — множество всех перестановок из т элементов; (—1)а — четность перестановки а е Бт.

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

Определение 3. Мультискобкой Кругликова-Лычагина

функций ^2,..., ^т+1 называется функция на пространстве джетов [2]

т+1

1, . . . , . . . , ¿т+1^ Fi,

г=1 ^

где Ьг = = (Ь1, Ь2,..., Ьт), г = 1,..., т+1, а ¿г означает пропуск ^й компоненты.

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

1) Если система Е формально интегрируема, то мультискобка

. . . , ^т+1}

равна нулю в силу системы Е.

2) Если система Е находится в полном пересечении ([1, 3]), то Е формально интегрируема тогда и только тогда, когда мультискобка

. . . , ^т+1}

равна нулю в силу системы Е.

Если рассматриваемая система состоит из более чем m + 1 уравнений, то необходимо рассмотреть скобки для всевозможных наборов из m+1 уравнений. Подробное изложение этих результатов можно найти в [5].

1. Вспомогательные процедуры

Для реализации вычисления формулы (5) нам необходимы следуюшие процедуры:

1) операции над дифференциальными операторами: сложение, умножение на функцию и композиция;

2) вычисление некоммутативного определителя матрицы;

3) вычисление линеаризации системы дифференциальных уравнений.

1.1. Операции над дифференциальными операторами

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

Таким образом, каждая пара ключ-значение ((cti, ..., an) ,Fa (x\,..., xn)) представляет оператор

Fa (xi,...,xn) d|a|/dxCT. Тогда весь ассоциативный массив T представляет дифференциальный оператор

где К (Т) — множество ключей массива. Далее в зависимости от контекста мы будем понимать под Т и сам оператор, и ассоциативный массив, которым оператор представляется.

aeK(T)

Сложение дифференциальных операторов. В выбранном представлении операторов сумма двух операторов

v-^ д|ст| v-^ д|т|

Tl = Е дХ^ и T2 = Е Gt дхт

CTGK(Ti) T€K(T2)

задается следующим образом:

дw

T+T2 = E ^+

aeK(Ti)\K(T2)

^ д|т | ^ д|ст| + E G ^ + E +G ) ^.

т €K(T2)\K(Ti) a€K(Ti)nK(T2)

Формула для разности аналогична. Реализация сложения двух операторов на языке Maple представлена ниже.

diffopPlus := proc( diffop1::table, diffop2::table )::table; local result::table, sigma::list; result := copy( diffop1 ); for sigma in indices( diffop2 ) do if sigma in {indices( result )} then

result[op( sigma )] := result[op( sigma )] + diffop2[op( sigma )]; else

result[op( sigma )] := diffop2[op( sigma )]; end if; end do;

return copy( result ); end proc;

Умножение дифференциального оператора на функцию.

Умножение оператора

Т = V

аЕК(Т)

на функцию д(жь ..., жп) слева задается формулой

д

E gF

aeK(T )

Код, ее реализущий:

diffopMult := proc( diffOp::table, g )::table; local result::table, sigma::list; result := copy( diffOp ); for sigma in indices( diffOp ) do

result[op( sigma )] := result[op( sigma )] * g; end do;

return copy( result ): end proc:

Композиция дифференцирования и дифференциального монома. Пусть дан моном

где т — мультииндекс. Рассмотрим композицию дифференцирования по переменной х% и монома М. Формула для вычисления композиции вытекает из правила Лейбница и выглядит так:

Процедура на Maple, вычисляющая данную композицию:

diffopDMonomial := proc( diffMon::table, varIndex::integer )::table; local result::table, sigma::list, othersigma::list; if nops( [indices( diffMon )] ) <> 1 then

error "Incorrect monomial has size %1, indices=%2", nops( [indices( diffMon )] ), [indices( diffMon )] end if;

if ivarsCount < varIndex or varIndex < 1 then

error "Incorrect variable's index %1", varIndex; end if;

sigma := indices( diffMon )[1]; othersigma := sigma;

othersigma[varIndex] := othersigma[varIndex] + 1; result := copy( diffMon );

result[op( sigma )] := TDiff( diffMon[op( sigma )], varIndex ); result[op( othersigma )] := diffMon[op( sigma )]; return copy( result ); end proc:

Композиция дифференцирования и дифференциального оператора. Используя функцию й1//орБМопот1а1, реализованную выше, мы можем построить композицию оператора взятия частной производной и дифференциального оператора. Пусть задан оператор

1

2

n

Ma,

a&C(T )

где Ma = F J. Формула для композиции записывается следу-

ющим способом:

Исходный код этой процедуры на Maple:

diffopDOperator := proc( diffOp::table, varIndex::integer )::table; local result::table, sigma::list, diffMon::table; if ivarsCount < varIndex or varIndex < 1 then

error "Incorrect variable's index %1", varIndex; end if;

result := table();

for sigma in indices( diffOp ) do

diffMon := table( [op( sigma ) = diffOp[op( sigma )]] ); result := diffopPlus( result, diffopDMonomial( diffMon, varIndex ) ); end do;

return result; end proc:

Композиция дифференциального монома и дифференциального оператора. Пусть

а1'1

N = 'aF

— дифференциальный моном, а

T = ^ Ma

a&C(T)

— дифференциальный оператор. Композиция N о T строится сде-дующим образом:

д|т I

N о T = N о Y, M = E N о M = f £ — о Ma.

ctGK(T) aEK(T) aEK(T)

Реализация на Maple:

diffopComposeWithMonomial := proc( diffMon::table, diffOp::table )::table; local result::table, diffs::list, f, diffIndex::integer, varIndex::integer; if nops( [indices(diffMon)] ) <> 1 then

error "Incorrect monomial has size %1, indices=%2", nops( [indices(diffMon)] ), [indices(diffMon)] end if;

result := copy( diffOp ); diffs := indices( diffMon )[1]; f := diffMon[op( diffs )]; for varIndex from 1 to ivarsCount do

for diffIndex from 1 to diffs[varIndex] do

result := diffopDOperator( result, varIndex ) end do end do;

result := diffopMult( result, f ); return copy( result ); end proc:

Композиция дифференциальных операторов. Имея в распоряжении композицию монома и оператора, мы можем легко построить композицию двух операторов. Пусть даны операторы

Т = ^ Ма

аеК(Тг)

и T2. Композицию Ti о T2 мы запишем в виде:

Ti о T2 = I Е ыа | о T2 = Е (M- о T2). VgK(Ti) / a€K(Ti)

Процедура, реализующая композицию двух операторов, представлена ниже.

diffopCompose := proc( diffOp1::table, diffOp2::table )::table; local result::table, sigma::list, diffMon::table; result := table();

for sigma in indices( diffOp1 ) do

diffMon := table( [op( sigma ) = diffOp1[ op(sigma) ]] ); result := diffopPlus( result,

diffopComposeWithMonomial( diffMon, diffOp2 ) ); end do;

diffopRemoveZeroes( result ); return copy( result ); end proc:

Применение дифференциального оператора к уравнению.

При реализации мультискобки нам необходимо применять операторы к уравнениям из исходной системы (1). Поэтому мы реализуем следующую процедуру.

Определим действие дифференциального оператора

T = V F д!

aeK(T)

на функцию G(xi, x2,..., xn) следующей формулой:

-(G)= Е ^

аек,(т)

Исходный код процедуры на Maple:

diffopOnEquation := proc( diffOp::table, eq ) local result, sigma::list, F, varIndex::integer, i::integer; result := 0;

for sigma in indices( diffOp ) do F := eq;

result := result + diffOp[op( sigma )] * TotDiff( F, sigma2vars( sigma ) ); end do;

return result; end proc;

1.2. Реализация некоммутативного определителя

Имея в распоряжении построенные вспомогательные процедуры, мы можем реализовать процедуру вычисления некоммутативного определителя в соответствии с определением 2. Мы будем вычислять определитель, раскладывая его по первому столбцу. Ниже мы приводим текст процедуры Ndet на языке Мар1е, вычисляющей определитель матрицы, используя рекурсивный алгоритм.

Ndet := proc( A::Matrix )::table;

local N::integer, temp::Matrix, result::table, row::integer, i::integer, j::integer, tempRow::integer, sign::boolean; N := op( 1, A )[1]; if N <> op( 1, A )[2] then

error "Matrix in Ndet is not square, %1 <> %2", N, op( 1, A )[2]; end if; if N = 2 then

return copy( diffopMinus( diffopCompose( A[1, 1], A[2, 2] ), diffopCompose( A[2, 1], A[1, 2] ) ) ); end if;

temp := Matrix( N - 1, N - 1 ); result := table(); sign := true; for row from 1 to N do sign := not sign; tempRow := 0; for i from 1 to N do if i <> row then

tempRow := tempRow + 1; for j from 2 to N do

temp[tempRow, j - 1] := A[i, j]; end do; end if; end do; if sign then

result := diffopMinus( result, diffopCompose( A[row, 1], Ndet( temp ) ) );

else

result := diffopPlus( result, diffopCompose( A[row, 1], Ndet( temp ) ) ); end if; end do;

return copy( result ); end proc:

Входной параметр функции Ndet имеет тип Matrix, который определен в Maple. Этот тип позволяет хранить матрицу, элементы которой могут иметь любой тип. В качестве результата работы функция возвращает дифференциальный оператор, хранящийся в ассоциативном массиве типа table.

1.3. Линеаризация

Важным этапом алгоритма вычисления мультискобки Кру-гликова-Лычагина является вычисление линеаризации системы

уравнений (1). Чтобы линеаризовать систему дифференциальных уравнений, необходимо линеаризовать каждое из уравнений этой системы. Используя явную формулу из определения (1), реализуем на языке Maple линеаризацию одного уравненения.

linearizeEquation := proc( F )::Array; local result::Array, i::integer, sigmas::set, sigma; result := Array( 1 .. dvarsCount ); for i from 1 to dvarsCount do result[i] := table(); sigmas := allSigmas( F, dvars[i] ); for sigma in sigmas do

result[i][op( sigma )] := diff( F, dvars[i][op( sigma )] ); end do; end do;

return result; end proc:

Аргумент функции linearizeEquation — это выражение F, содержащее независимые переменные (хранятся в глобальном списке ivars) и координаты джетов, записанные через зависимые переменные (хранятся в глобальном списке dvars) с мультииндек-сом. Эта функция возвращает объект типа Array (определенный в Maple тип массива), содержащий дифференциальные операторы.

Размер возвращаемого массива в точности равен числу m — количеству неизвестных функций Ui, которое равно переменной dvarsCount — длине списка dvars. На i-й позиции в этом массиве стоит дифференциальный оператор при i-й функции в формуле (3).

Встречающася в приведенном тексте функция allSigmas — одна из вспомогательных функций, на реализации которых мы не останавливаемся. Данная функция, анализируя выражение, передаваемое ей в первом аргументе F, ищет все записанные джета-ми производные функции, имя которой было передано ей вторым аргументом dvars[i], и возвращает список всех мультииндексов, стоящих при этом имени в выражении F.

Построенную линеаризацию уравнения мы применим для процедуры линеаризация системы уравнений вида (1):

linearizeSystem := proc( sys::list )::Matrix;

local result::Matrix, eq::Array, i::integer, j::integer, eqsCount::integer; eqsCount := nops( sys );

result := Matrix( eqsCount, dvarsCount ); for i from 1 to eqsCount do

eq := linearizeEquation( sys[i] );

for j from 1 to dvarsCount do

result[i, j] := copy( eq[j] ); end do; end do;

return result; end proc;

Процедура linearizeSystem линеаризует систему уравнений, переданную ей в списке sys. Находя линеаризацию каждого уравнения в этом списке, она строит матрицу дифференциальных операторов, имеющую dvarsCount столбцов и на единицу большее, чем число уравнений, количество строк.

2. Реализация мультискобки Кругликова-Лычагина

Реализовав все вспомогательные процедуры, мы переходим к реализации определения (3). Ниже представлен исходный текст процедуры MultiBracket.

MultiBracket := proc( eqSys::list )

local linearization::Matrix, Ndets::Array, result, i::integer; linearization := linearizeSystem( eqSys ); Ndets := allNdets( linearization ); result := 0;

for i from 1 to dvarsCount + 1 do

result := result + (-1)Л( (i + 1) mod 2 ) * diffopOnEquation( Ndets[i], eqSys[i] ); end do;

return simplify( result ); end proc;

Входной параметр eqSys — список уравнений исходной системы. Процедура возвращает выражение — вычисленную муль-тискобку. Процедура allNdets принимает на вход матрицу дифференциальных операторов размера m х (m + 1), а возвращает список из m некоммутативных определитей, полученных из входной матрицы.

Замечание 1. Оценим сложность вычисления нескольких реализованных функций. Функция Ndet вычисляет определитель матрицы разложением по столбцам, следовательно, выполняет m! операций композиции, где m — число неизвестных функций. Функция MultiBracket вызывает m + 1 раз функцию Ndet, а затем применяет результат каждого вычисления m + 1 раз уравнениям исходной системы, т.е. выполняет (m + 1)2 m! композиций.

Несмотря на то, что формула вычисления мультискобки выполняет значительное количество операций, этот метод получения условий формальной интегрируемости является довольно быстрым, потому что задает явную формулу. Сравнение метода скобок с другим методом - дифференциальными базисами Гребнера -приведено в [4].

3. Пример

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

х = а(х,у) дх+Ъ(х,у) ду.

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

{ Ап = 0,

| X(п) = 0,

где А = дхх + дуу — оператор Лапласа. Согласно условию Коши-Римана эта система эквивалентна следующей системе дифференциальных уравнений первого порядка:

{Пх - Уу = 0,

Пу + Ух = 0,

а(х,у)пх + Ъ(х,у)пу = 0. В силу того, что векторное поле X не имеет особых точек, можно считать, что функция а(х, у) нигде не обращается в ноль. Тогда без ограничения общности можно положить а(х, у) = 1. Система

(6) примет вид:

Пх - Уу = 0,

(7) { Пу + Ух = 0,

Пх + Ъ(х,у)Пу = 0.

Программа на Maple, вычисляющая скобку для системы (7):

# подключаем модуль Brackets with( Brackets ):

# настраиваем внутренние параметры этого модуля:

# списки независимых и зависимых переменных setup( [x, y], [u, v] ):

# задаем систему уравнений

S :- [ u[1, 0] - v[0, 1], u[0, 1] + v[1, 0], u[1, 0] + b(x, y) * u[0, 1] ]:

# вычисляем скобку и группируем члены у проивзодных u_x, u_y, u_xy MlBr :- collect( MultiBracket( S ), [ u[0, 1], u[1, 1], u[0, 2] ] );

В результате компьютерных вычислений получим мультис-кобку Кругликова-Лычагина:

(8) М = АЬиу + 2Ьу иуу + ,

где А = + — оператор Лапласа. Ограничим скобку М на систему (7), используя третье уравнение = —Ь(ж, . Дифференцируя его по у, получим = —— . Ограниченную скобку обозначим

М£ = (АЬ — 2ЬжЬу )■% + 2(Ьу — .

Чтобы система (7) была формально интегрируема, необходимо и достаточно, чтобы М^ = 0. Видно, что это условие выполнено тогда и только тогда, когда коэффициенты при и равны нулю. Это условие приводит к следующей системе уравнений:

( АЬ = 2Ьу

(9) Ь ' "

Решая ее, получим

Ь(ж,у) = Л-1 (у)ж + ^2(у), где ^1, Л,2 — произвольные гладкие функции.

Чтобы исследовать совместность системы (7), перепишем ее в форме, не содержащей функцию -и(ж, у), и добавим к ней уравнение Ме = 0:

+ — 0

(10) <( + = 0, (АЬ — 2ЬхЬу )иу + 2(Ьу — = 0.

Продифференцируем второе уравнение этой системы по х и по у

и добавим полученные соотношения к системе (10):

+ - 0

их + Ьиу — 0,

ижж + Ьжиу + Ьижу - 0,

ижу + Ьу Иу + Ьиуу — 0,

^ (АЬ — 2ЬхЬу)иу + 2(Ьу — ЬЬх)и.

уу

0.

Эта система эквивалентна следующей системе:

ижж + иуу — 0,

+ Ьиу — 0,

иху + Ьу Иу + Ьиуу — 0, — (1 + Ь2)иуу + (Ьх — ЬЬу )иу — 0,

0.

(11)

ч 2(Ьу — ЬЬх)иуу + (АЬ — 2ЬхЬу)иу Необходимо рассмотреть два случая.

В первом случае, когда иу — 0 и иуу — 0, определитель системы, составленной из четвертого и пятого уравнений системы (11), должен быть равен нулю. Это условие дает следующее соотношение на функцию Ь(х, у): (12) (Ь2 + 1)АЬ — 2Ь(ЬХ + Ь2) — 0.

C учетом соотношения (12) выразим из системы (11) производные ижж и иху через иу:

(Ьх — ЬЬу) Иу

—

Ь2 + 1

и

ху

— — Ьу +

Ь (Ьх — ЬЬу) Ь2 + 1

Иу.

Используя эти соотношения, условие (12) и систему (11), можно показать, что условия формальной интегрируемости системы (11):

И

жжу

И

жуж

выполняются тождественно.

иуух - ижуу

Рассмотрим второй случай. Пусть = = 0. Тогда система (11)сводится к системе:

i = 0

\ Ux = 0,

решение которой u(x, у) = const. Таким образом мы приходим к следующей теореме.

Теорема 1. Векторное поле

^ д w ч д

X = дХ + b(x,y) ду

обладает не равным константе гармоническим первым интегралом тогда и только тогда, когда функция b(x, у) удовлетворяет следующему дифференциальному уравнению: (b2 + 1)Ab - 2b(bX + ) = 0.

Литература

1. АЛЕКСЕЕВСКИЙ Д.В., ВИНОГРАДОВ А.М., ЛЫЧА-ГИН В.В. Основные идеи и понятия дифференциальной геометрии // Современные проблемы математики. Фундаментальные направления. - 1988. - Т. 28 (Итоги науки и техники ВИНИТИ АН СССР). - С. 5-289.

2. ВИНОГРАДОВ А.М., КРАСИЛЬЩИК И.С., ЛЫЧА-ГИН В.В. Введение в геометрию нелинейных дифференциальных уравнений. - М.: Наука, 1986. - С. 336.

3. КОКС Д., ЛИТТЛ ДЖ., О'ШИ ДЖ. Идеалы, многообразия и алгоритмы.. — М.: Мир, 2000. — С. 687.

4. KRUGLIKOV B. Note on two compatibilty criteria: Jacobi-Mayer bracket vs. differential Groebner basis // Lobachevskii Journal of Mathematics. - 2006. - Vol. 23. - С. 57-70.

5. KRUGLIKOV B., LYCHAGIN V. Multi-brackets of differential operators and compatibility of PDE systems // C. R. Acad. Sci. Paris. - 2006. - Vol. 342. - С. 557-561.

6. Maple Advanced Programming Guide / Maplesoft, a division of Waterloo Maple Inc. 1996-2009. -

URL: http://www.maplesoft.com/view.aspx?sl=32470 (дата обращения: 12.02.2011).

IMPLEMENTATION OF THE KRUGLIKOV-LYCHAGIN MULTIBRACKET ON MAPLE

Sergey Tychkov, Institute of Control Sciences of RAS, Moscow, postgraduate student, ([email protected]).

Abstract: We describe implementation of Kruglikov-Lychagin multibracket on Maple. An example is given of application of Kruglikov-Lychagin multibracket for compatibility analysis of a specific system of partial differential equations.

Keywords: differential equations, Kruglikov-Lychagin multibracket, Maple.

Статья представлена к публикации членом редакционной коллегии В. А. Уткиным