CC BY
70
УДК 517.55, 519.6
М. В. Бураченко
АЛГОРИТМ ПОНИЖЕНИЯ КРАТНОСТИ ИНТЕГРАЛА РАЦИОНАЛЬНОЙ АЛГЕБРАИЧЕСКИ ТОЧНОЙ ДИФФЕРЕНЦИАЛЬНОЙ ФОРМЫ С КВАЗИЭЛЛИПТИЧЕСКИМ ЗНАМЕНАТЕЛЕМ
Изложен метод понижения кратности интеграла рациональной алгебраически точной дифференциальной формы с квазиэллиптическим знаменателем вписан алгоритм и представлена его реализация для системы компьютерной алгебры Maple 9.
В данной статье в качестве основного объекта исследования будем рассматривать интеграл вида
I =| Р(х) Сх =| Р(Хп) dx1...dxn , (1)
RnQ( х) RnQ( Х1,..., Хп) 1 п
где Р(х) и Q(x) - полиномы. Предполагается, что Q не обращается в нуль на Яп и интеграл (1) сходится абсолютно. Такие интегралы встречаются в различных областях математики и теоретической физики. В частности, к ним приводятся фейма-новские интегралы в квантовой теории поля [1], [2]; некоторые интегралы такого типа выражают топологические заряды инстантонных полей [3].
В работе [4] была предложена формула для понижения кратности подобных интегралов в случае, когда Q(x) - эллиптический полином, не обращающийся в нуль на Яп, а подынтегральная форма алгебраически точна. В работе [5] этот результат был обобщен на случай квазиэллипти-ческого знаменателя.
В данной статье, развивающей идеи работы [5], излагается метод понижения кратности интеграла (1) в случае квазиэллиптического знаменателя, описывается соответствующий алгоритм и представляется его реализация для системы компьютерной алгебры Мар1е 9, позволяющей вычислять интеграл (1) до конца в пространстве Я2 и понижать кратность интеграла на единицу в пространстве Я3.
Основные определения и постановка задачи
Введем, следуя [5], некоторые определения и условные обозначения, необходимые для формулировки основного результата.
Ниже будут рассмотрены интегралы вида (1), где Р и Q - полиномы с рациональными коэффициентами. Это ограничение на коэффициенты связано с возможностью реализации алгоритмов. С теоретический точки зрения допустимы коэффициенты из любого упорядоченного поля характеристики нуль.
Пусть полином Q(x) представлен в виде
^Х)= Т Са Х“ = Т Са1...ап Х1 “'...'О ,
а<^Ып (,...,ап )ып
где N - множество натуральных чисел с добавленным нулем. Носителем Q будем называть множество suppQ = {а е N | са Ф 0}.
Многогранником Ньютона Д = Д^), соответствующим полиному Q, называется выпуклая оболочка в Яп для носителя Q.
Срезкой Qv полинома Q в направлении ковек-тора V е Яп (здесь Яп означает сопряженное пространство к Яп) будем называть полином
Qv = Тсаха , где Д” = {^ еД: (м>, V^ = тт(м,V)} -
аеДV МеД
грань Д в направлении V.
Определение 1. Полином Q называется квазиэллиптическим, если для любого ненулевого ко-вектора V є Я"* срезка Qv не обращается в нуль в торе (Я \ {0})".
Замечание 1. Так как число граней многогранника Д конечно, то условие квазиэллиптичности достаточно проверять лишь для конечного числа срезок Qv.
При дальнейшем изложении будем полагать, что для многогранника Ньютона Д, соответствующего полиному Q, справедливо следующее предложение:
Предположение 1. Многогранник Д(О) содержит по ребру на каждой координатной оси Ох,, 1 < і < " .
Каждому многограннику Ньютона А можно поставить в соответствие двойственный конический полиэдр (веер) ХД. Дадим его определение.
С каждым ковектором V є Я"* связана грань Д1” многогранника Д, на которой скалярное произведение ^, и> минимально на Д. Два ковектора будем считать эквивалентными, если связанные с ними грани совпадают. Замыкание классов эквивалентности ковекторов, связанных с гранью Д,, образует конус в Я"*, порожденный элементами решетки . Этот конус называется двойственным к грани Ді. Совокупность конусов, двойственных всем граням многогранника Д, образует веер ХД, который называется двойственным к многограннику Д. Мы будем рассматривать только полные вееры, для которых объединение вхо-
Г)"*
дящих в них конусов есть все пространство Я .
Очевидно, что образующими ребер веера ХД являются внутренние нормали Vj гиперграней
А("-1) є А , а набор V71,..., v]k таких нормалей порождает конус ст(у71,..., Vі7) в ХД тогда и только
тогда, когда соответствующие грани
Нхи
Д("-1) Д("-1)
А 71 ,..., А 7к
имеют непустое пересечение
А71 7^. Размерность конуса о будет дополнительной к размерности грани А^ , т. е.
= п-^шД. . .
Конус, число ребер которого равно его размерности, назовем симплициальным. Симплици-
альный к-мерный конус ст(,...,Vі)
примитивным, если
называется
параллелепипед
П=і%XvV1,0<X< 11 не содержит целых точек
(точек, все координаты которых целые числа), отличных от точки 0, или, что то же самое, векто-
1 к /Т
ры V , .., V можно дополнить до базиса решетки . При к = п это условие эквивалентно тому, что
гг 1 к
матрица, образованная векторами V , .., V - уни-модулярна, т. е. ее определитель равен ±1. Веер будем называть примитивным, если любой его конус симплициален и примитивен.
Пусть ё - число гиперграней многогранника Д^) и ХД - веер, двойственный А^) . Обозна-1 ё
чим через V ,...,V минимально целочисленные образующие ребер веера ХД. Если веер ХД не является примитивным, то можно построить примитивный веер ХА , состоящий из примитивных конусов, причем каждый конус ХА лежит внутри некоторого конуса из Ед , и если конус из ХД - примитивный, то он не измельчается. Пусть в резуль-
гг 1 ё
тате такого подразбиения к векторам V , ..., V дог ё+1 ё+г
бавились векторы V , ..., V .
Для каждого вектора V , 1 < к < ё, определим
матрицу Ск = Ск1 .кп-1к = (vkl,..., ^"-1,^), столб-
__ к кп 1
цами которой служат векторы V 1,..., V п 1, порождающие вместе с Vі конус в ХА. Поскольку ХА - примитивный веер, то все матрицы Ск унимо-дулярны, т. е. det Ск = ±1, 1 < к < ё.
С каждой матрицей Ск свяжем биномиальную Ск
замену координат х = у к, что в терминах координат (^к,..., v"<:) векторов Vі
V" I векторов V означает следующее:
V, ^ v1
х1 = У11 ... У"-1 У"
^1 ^И-1 vк
х" = У1Я ...У"-1 Уп .
На подынтегральное выражение в (1) будем смотреть как на дифференциальную форму, и записывать ее в виде
> х )
—’Цгёхх л... л ёхп
(2)
Q{Х1, ..., Хп)
Определение 2. Дифференциальная форма ю называется алгебраически точной, если существует рациональная (п - 1) - форма ф, дифференциал Сф которой равен ю. Такую форму ф будем называть рациональной первообразной для формы ю.
Определение 3. Пусть ф - рациональная форма степени ^, т. е. ф= Т Ф,,..,р (х)Сх,1 л... л Сх,г ,
1</'1<...</р <п
где коэффициенты фг[ .р (х) - рациональные
функции переменной х е Яп . Если знаменатели всех функций ф. . являются степенями квазиэл-
т 1\...1р
липтического полинома Q, то мы говорим, что форма ф не растет на бесконечности, когда для каждой гиперграни Д— многогранника Д^) форма ф(уС )
не имеет полюса на гиперплоскости уп = 0.
Введем
обозначения:
V = V1 + ... + V"
0± =](,... У"-1 )є Я"
: У1
Vі"-1 -1 >
• у"-1 0^.
В указанных обозначениях справедлива следующая теорема:
Теорема 1. Пусть Q - квазиэллиптический полином, не обращающийся в нуль на Я", и интеграл (1) абсолютно сходится. Если дифференциальная форма (2) имеет рациональную первообразную ф, которая не растет на бесконечности, то для интеграла (1) справедлива формула
1 = (-1)п% + (-1)Ы)^Ск'|[ф(УСк к=о, (3)
к=1 ^ о+ о_ )
где суммирование ведется по всем (п - 1)-мерным граням д"- (Q).
Доказательство этой теоремы в рамках данной работы не приводится, но его можно найти в [5].
Для случая п = 2 формулу (3) можно переписать в виде
1=2, % ЛрИ )]іу
к: V -четн Я
|ІУ2 =0>
Алгоритм понижения кратности интеграла рациональной алгебраически точной дифференциальной формы с квазиэлептическим знаменателем
На основании вышеизложенного представим алгоритм понижения кратности интеграла (1). Все вспомогательные алгоритмы, выделенные латинскими буквами, будут описаны далее, после основного алгоритма.
I=1
<
QuasiellipticIntegration(P, Q, x)
Input: P = P(xi, ... xn) - полином, зависящий от n переменных, - числитель подынтегральной формы интеграла (1);
Q = Q(x1, ... xn) - полином, зависящий от n переменных, - знаменатель подынтегральной формы интеграла (1);
x = [x1,.xn] - список переменных;
Output: интеграл вида (3), кратность которого (n - 1);
1) проверка того, что Q не обращается в
Г)П
нуль в R :
t1 : = PolyNot0(Q, x);
if t1 = false then error 'знаменатель обращается в нуль';
2) проверка квазиэллиптичности полинома Q:
t2 : = Quasielliptic(Q, x);
if t2 = false then error 'знаменатель не квазиэллиптический';
3) проверка условия предположения 1:
Д : = NewtonPolyhedron(Q, x);
if найдется 1 < I < n, такой что не существует 5 е N, такого что I о,..., s, ...,0|еЛ then error
с. = (v vkl
vk"--l. vk
условие предположения 1 не выполняется';
4) проверка того, что дифференциальная P
форма с коэффициентом q алгебраически точная, и вычисление ее рациональной первообразной ф, не растущей на бесконечности:
ф : = RationalPrimitive (P, Q, x) ф: = = RationalPrimitive (Q, P, x);
if ф = [О,...,О] then error 'для исходной
n
подынтегральной формы не существует рациональной первообразной, не растущей на бесконечности';
З) проверка абсолютной сходимости интеграла (1):
t3 = AbsoluteConvergence(P, Q, x); if t3 = false then error 'интеграл не сходится абсолютно';
б) Х : = NormalFan(Q, x) - веер в R"*, двойственный многограннику Д, конуса которого порождены минимальными целочисленными ко-векторами vk, k = l,...,d внутренних нормалей к (n - 1)-мерным граням A"-1 многогранника Д;
У) ~ : = PrimitiveNormalFan(X) - примитивный веер, полученный в результате простого подразбиения веера Х (Q) добавлением к его образующим ковекторов vd+1 кvd+r;
В) каждому ковектору v , k = 1, ..., d, сопоставим унимодулярную матрицу
), так чтобы ковекторы v",...,vkn-1,vk образовывали конус в X(Q);
9) с каждой матрицей Ck, 1 < k < d + r свяжем биномиальную замену переменных
x = yCk;
10) return (интеграл кратности (n - 1), записанный по формуле (2))
end.
Пока интеграл, получившийся в результате работы алгоритма, удовлетворяет условиям теоремы 1 и предположения 1, его кратность можно понижать, используя тот же алгоритм. Но если итоговый интеграл не будет удовлетворять какому-либо из условий, то его можно попытаться вычислить стандартными средствами.
Разберем подробнее вспомогательные алгоритмы, встречающиеся в основном алгоритме, и особенности их реализации в системе компьютерной алгебры Maple.
Построение многогранника Ньютона и двойственного ему веера. При рассмотрении Quasiel-lipticIntegration, в частности при проверке условий квазиэллиптичности знаменателя, алгебраической точности подынтегральной формы и абсолютной сходимости интеграла, нам понадобятся алгоритмы построения многогранника Ньютона (алгоритм NewtonPolyhedron) и двойственного ему веера нормалей (алгоритм NormalFan). Напомним, что многогранником Ньютона полинома называется выпуклая оболочка точек носителя этого полинома.
Алгоритмы построения многогранника Ньютона полинома и двойственного ему веера осуществляются в два этапа: на первом этапе определяется носитель полинома, а на втором находится его выпуклая оболочка.
Для построения выпуклой оболочки множества точек n-мерного пространства используется алгоритм, базирующийся на алгоритме Моцкина-Бургера [6], и его реализация П. А. Буровским [7] в виде Maple-программ Convexhull и FullConvex-hull программного пакета convexhull.
Программа Convexhull по заданному списку точек и размерности пространства строит список точек исходного множества, попадающих на каждую из гиперграней выпуклой оболочки точек исходного списка, и список образующих нормальных конусов соответствующих гиперграней. Программа FullConvexhull строит описание граней всех размерностей и нормальных конусов к ним.
Для построения многогранника Ньютона и двойственного ему вектора нормалей (алгоритмы NewtonPolyhedron и NormalFan соответственно)
Зб
мы будем пользоваться какой-либо одной из этих процедур в зависимости от полноты требуемых данных.
Если исходный интеграл берется по пространству R2, то подключение пакета convexhull не требуется. Для построения многогранника Ньютона (это будет выпуклый многоугольник) используется стандартная Maple-процедура convexhull программного модуля simplex. Для построения веера нормалей векторы, получающиеся при обходе многоугольника против часовой стрелки, поворачиваются на угол 900 (если (x, у) - координаты вектора-стороны многоугольника, то (-y, x) - координаты внутренней нормали к этой стороне) и минимизируются т. е. для каждого вектора вычисляется наименьший общий делитель его координат и все координаты на него делятся.
Проверка того, что знаменатель не обращается в нуль, и квазиэллиптичность знаменателя. На знаменатель Q(x) подынтегральной формы интеграла (l) накладываются два условия:
Условие 1. Полином Q(x) не обращается в
П"
нуль в R .
Условие 2. Полином Q(x) является квазиэллип-тическим полиномом, т. е. для любого ненулевого ковектора v є Rn* срезка Qv не обращается в нуль в торе (R \ {0})".
Для проверки этих условий используется так называемый алгоритм цилиндрического алгебраического разбиения, или сокращенно CAD-алгоритм (от англ. Cylindrical Algebraic Decomposition). Этот алгоритм является инструментом для решения систем полиномиальных урав-
^ п"
нений и неравенств в R и заключается в построении такого разбиения пространства R", чтобы в каждой компоненте разбиения сохранялись знаки всех полиномов исходной системы. Под разбиением пространства будем понимать конечную совокупность непересекающихся областей пространства, объединение которых равно этому пространству
Впервые такой метод был предложен в 194В г. Тарским, но он не был реализован практически. В 19У3 г. Коллинз предложил конструктивный метод. В рамках данной статьи алгоритм был реализован в виде Maple-программы CAD.
Введем некоторые определения, необходимые для формулирования основных положений метода. Все определения даются согласно работам [В, 9]. Более детальное описание алгоритма можно найти в [9].
Входными данными для алгоритма является множество полиномов Фс Q X.,..., X" ], где
Q - поле алгебраических чисел, и список переменных [x.,..., xn ].
CAD-алгоритм состоит из трех основных этапов: проекции, построения базы и расширения.
Этап проекции состоит из (n - 1)-го шага. На каждом шаге строится новое множество полиномов, зависящих от числа переменных, на единицу меньшего, чем полиномы множества, полученного на предыдущем шаге. Множества нулей полиномов, полученных на каждом шаге этапа, - это проекции особых точек множеств полиномов предыдущего шага, т. е. точек самопересечений, изолированных точек и точек, в которых касательные вертикальны. Переход от множества полиномов к множеству их проекций осуществляется при помощи оператора proj. Чтобы его записать, нам понадобятся некоторые условные обозначения.
Пусть F = {/1, ...,f}, гдеf(xb ..., x„)eOh,...,xn], 1 < i < r. Эти полиномы можно рассматривать как полиномы, зависящие от одной переменной xn, коэффициентами которых являются полиномы от (n - 1)-й переменной x1, ..., xn-1, т. е. fi(x1, ..., xn), 1 < i < r, представляется в виде f(x1,..., xn-1,xn) =
= fi,di Xx1,..., xn-1Yn + ... + f ,0 Xx1,..., xn-1 ) е Q Х — xn-1 k ]
и fik (x1,.., xn-1) - коэффициент i-го полинома при
k-й степени переменной, xn - полином, зависящий от переменных x1, ., xn-1.
Обозначим fk (, ..., xn )= £fj (*1,..., xn-1 ) ,
j=0
где 0 < k < dj.
Обозначим символом Dxоператор дифференцирования по переменной xn.
Матрицу mat(gb ., gk) = [g i-g] размера k x k ,
где l = 1 + max (deg(gi)), назовем матрицей, ассо-
1<i<k
циированной с полиномами g1, ., gk.
Для полинома g степени s по переменной xn и полинома h степени m по переменной xn определим psc (g, h)= det*Ml), где матрица Mi состоит из l первых столбцов матрицы mat (xs-l-1 g,..., g,xm-l-1h,..., h), ассоциированной с
полиномами xs-1 -1 g,..., g,xm-1 -1h,..., h .
Тогда
proj(F) = proj^F) и proj2(F) и proj3(F), где proj1 = {fik(x1,..., xn-1) 11 < i < r, 0 < k < dt};
proj2 = {psc Xln *^k (x1,..., xn), Dxn (x1,..., xn)))
|1 < i < r, 0 < l < k < di};
proj3 = {psc m i,fIki (xl,..., xn), fk}j (xl,..., xn))
|1 < i < j < r, 0 < ki, kj < dt ,0 < m < min(ki, kj ) < di}.
Определенный таким образом оператор proj рекурсивно применяется (n - 1) раз к исход-
ЗУ
ному множеству полиномов, т. е. рщ)0(Ф) = Ф, рго^Ф) = рго)(рго)0(Ф)), ..., ргоУ(Ф) = ргоіСршТ^Ф)), ..., proj"-1(Ф) = ргсу(ргсу"-2(Ф)). Последнее множество - множество полиномов, зависящих от одной переменной х1, - конечный результат этапа проекции.
На этапе построение базы строится множество пробных точек разбиения пространства Я1. В качестве входных данных выступает множество полиномов рщ)"-1(Г) - это множество полиномов, зависящих от одной переменной, получившееся в результате этапа проекции, и сама переменная, от которой зависят данные полиномы. Вычисляются действительные корни полиномов множества рщ)"-1(Г). Эти корни в общем случае являются алгебраическими числами и могут быть представлены в виде пары полинома и интервала, на котором этот полином имеет единственный корень -данное алгебраическое число, причем концами интервала являются рациональные числа, а полином выбирается таким образом, чтобы его степень была наименьшей. Пусть у нас получилось следующее множество алгебраических чисел ?1 є(и1,V!] к, 5, є (и,vs], где все и7, vj є ё. про6-ные точки выбираем следующим образом: а1 = и1, а2 = £ь..., а2х+1 = vs + 1. Последнее множество и является результатом этапа построения базы.
Цель этапа расширения заключается в построении разбиения пространства Я". В качестве входных данных выступает множество пробных точек, полученное на этапе построения базы: Ьа8е0 = {а1, ..., а2а+1}. Этап расширения содержит (п - 1) шагов.
На первом шаге точки множества Ьа8е0 под______________________________ • "-2 /т-\
ставляются в полиномы множества ргс) (г) - про-
екции, полученной на (п — 2)-м шаге этапа проекции. В результате такой подстановки для каждой точки а7, 1 < 7 < 2s + 1 получается множество полиномов, зависящих от одной переменной х2. К ним снова применяется этап построения базы. Пусть в результате для каждой точки а7, 1 < 7 < 2s + 1 получилось множество пробных точек 1,..., Р2г +1}. Тогда множество пробных точек
разбиения Я2 организуется следующим образом:
На втором шаге эти точки подставляются в полиномы из множества рщ)"-3(Г), и аналогичная процедура повторяется. Таким образом, после выполнения (п — 1)-го шага мы получим множество пробных точек пространства Я". Это и есть результат СЛБ-алгоритма.
Корректность рассмотренного СЛБ-алгоритма в рамках данной статьи не доказывается. Это доказательство можно найти, например, в [9].
Вернемся теперь к проверке условий 1 и 2. Достаточным условием того, что полином не обращается в нуль в Rn, является то, что все входящие в его состав мономы с xf1 к xa'n имеют по-
а 1 n
ложительные коэффициенты са и четные степени ai по каждой переменной xi, а свободный член полинома не равен нулю. Если это условие не выполняется, то для проверки условия 1 используется CAD-алгоритм. С его помощью строится множество пробных точек для этого полинома; затем вычисляются значения исходного полинома в этих точках, и если ни в одной из этих точек полином не обращается в нуль, то можно заключить, что полином не обращается в нуль в Rn. Если же хотя бы в одной пробной точке значение полинома равно нулю, то условие 1 не выполнено.
Таким образом, алгоритм PolyNot0 выглядит следующим образом:
PolyNot0(Q, x)
Input: Q = Q(x1, ..., xn) - полином, зависящий от n переменных;
x = [ x1, ., xn] - список переменных;
Output: true, если полином Q не обращается в нуль в Rn, false - в противном случае;
if (все мономы полинома Q имеют положительные коэффициенты и четные степени по каждой переменной) and (свободный член полинома Q равен нулю) then return (true) else
1) SamplePoints : = CAD({Q}, x) - множество пробных точек разбиения пространства Rn;
2) ValueQ - множество значений полинома Q в точках множества SamplePoints;
3) if 0 не содержится в множестве ValueQ then return (true) else return (false) end if
end if end.
Достаточным условием квазиэллиптичности полинома является то, что все входящие в его состав мономы с xf1 к xa'n имеют положительные
а 1 n
коэффициенты c и четные степени i по каждой переменной xi. Если это условие не выполняется, то при помощи CAD-алгоритма нужно построить множество пробных точек для всех срезок этого полинома (не имеет значения, рассматривать эти полиномы вместе или по отдельности), исключить из этого множества точки, у которых хотя бы одно координата равна нулю (это будет соответствовать исключению координатных осей из R , так как в данном случае нас интересует только тор (R \ {0})”), подставить оставшиеся точки в срезки исходного полинома и убедиться, что они не обращаются в нуль ни в одной точке из оставшихся точек. Алгоритм Quasielliptic выглядит таким образом:
ЗБ
Quasielliptic(Q, x)
Input: Q = Q(x., ..., X") - полином, зависящий от n переменных;
x = [ x., ., xn] - список переменных; Output: true, если полином Q квазиэллиптиче-ский, false - в противном случае;
if (все мономы полинома Q имеют положительные коэффициенты и четные степени
по каждой переменной) then return (true) else
1) A(Q):=NewtonPolyhedron(Q,x) - многогранник Ньютона полинома Q (здесь нас будут интересовать все грани многогранника Ньютона, поэтому мы будем использовать процедуру FullConvexhull);
2) F = {/і, ..., fk} - множество срезок, соответствующих граням многогранника A(Q);
3) SamplePoints : = CAD({f1, ... f},x) - множество пробных точек разбиения пространства R";
4) SamplePointsInTor - подмножество точек из SamplePoints, у которых ни одна координата не равна нулю;
З) ValueQ - множество значений полиномов множества F в точках множества Sample-Points;
б) if 0 не содержится в множестве ValueQ then return (true) else return (false) end if end if end.
Вычисление рациональной первообразной, не растущей на бесконечности. Алгоритм нахождения рациональной первообразной (Rational-Primitive), не растущей на бесконечности, основывается на методе неопределенных коэффициентов Остроградского и осуществляется в три этапа:
- определяются степени мономов, входящих в состав коэффициентов первообразной;
- выписываются полиномы с неопределенными коэффициентами и найденными степенями;
- решается уравнение на коэффициенты.
В качестве входных данных для этого алгоритма выступают полиномы P(x), Q(x) и список переменных x. Если исходная дифференциальная форма алгебраически точна и имеет не растущую на бесконечности рациональную первообразную, то алгоритм возвращает значение этой первообразной, если нет, то возвращается список из n нулей.
Рассмотрим теперь все этапы алгоритма более подробно.
Допустим, что дифференциальная форма ю = •••, Xj)) к,dxn алгебраически точна и
-, Xn ) "
ф - ее рациональная первообразная, не растущая на бесконечности. Представим ю в виде
dx следующим образом: если найдется
Q2k (x) = Q(x), то тогда P(x) = P(x); если найдется полином Q(x) и натуральное число k > 1, такие что Q2k (x) = Q(x), то тогда p(x) = p(x)^ Q(x); если таких Q(x) и k не существует, то тогда
q(x) = q(x) , k = . P(x) = P(X) • Q(x) .
Рациональную первообразную ф будем искать в виде
1 "
Ai (x)dxi л к [i]- к л dX" ,
B(x )i=i
где x = (x., ..., X"), A, - некоторые полиномы; B = Q2k-1 - квазиэллиптический полином, не об-
^ Г»"
ращающийся в нуль на R .
Пусть A = A(B) - многогранник Ньютона. Очевидно, что многогранники A(Q) и A(B) гомотетичны. Таким образом, если A(Q) удовлетворяет условию предположения l, т. е. содержит по ребру на каждой координатной оси в R", то и A(B) удовлетворяет этому условию. Веера нормалей, двойственные к A(Q) и A(B), совпадают. Пронумеруем векторы веера нормалей, двойственного
A(B) таким образом, чтобы vi =| 0, .„,1,к,0 |єR"*,
P(x) Q* (x)
полином q(x) и натуральное число k, такие что
і = 1, ..., ".
Обозначим через Д множество внутренних точек многогранника Д, а через А0к - множество точек, лежащих в относительной внутренности
(п -1)-мерной грани Ак = А*"-1^), двойственной
к
вектору V.
В работе [10] доказывается теорема о том, что в предположении 1 для того чтобы дифференциальная форма ф не росла на бесконечности, ее коэффициенты должны иметь вид
Д(х) = (-1)х,. ЬкДк(х), I = 1, ..., п,
к=п+1
где Ді0 и Дк - полиномы, многогранники Ньютона которых удовлетворяют условиям
а(д0)+1, с А0 и Ао, А(Дк)+1 сА0к , здесь
I = (1, к,1)єЯ", I. =^1, ^,0, к1^є Я", а Vlk , ..., V"
координаты вектора, двойственного грани
А"- (В ).
Исходя из этой теоремы, определим возможные степени полиномов Д,(х), I = 1, ..., ", коэффициенты которых не определены. Затем эти полиномы подставляем в уравнение
% 1)і-1 [В(х)дДхх - (2к - 1)Ді (х))ВМ| - Р(х) = 0
і=1 ^ дхі дхі )
и решим уравнение на коэффициенты.
Если в результате решения этого уравнения все неопределенные коэффициенты получились равными нулю, то рациональной первообразной ф, не растущей на бесконечности, не существует,
и мы возвращаем список из n нулей; если же не все неопределенные коэффициенты равны нулю, то возвращается список коэффициентов первообразной.
Абсолютная сходимость интеграла. В [2] доказывается следующая теорема:
Теорема 2. Если Q - квазиэллиптический полином, не обращающийся в нуль на R", то интеграл (1) абсолютно сходится тогда и только тогда,
когда I + A(p) с A0 (Q), т. е. когда сдвиг многогранника A(P) вдоль вектора I = (l, ...,і)єR" лежит во внутренности A0(Q) многогранника A(Q).
На ее основе составлен следующий алгоритм:
AbsoluteConvergence(Q, P, x)
Input: P = P(xi, ..., X") - полином, зависящий от n переменных, числитель подынтегральной формы интеграла (l);
Q = Q(x1,.,xn) - полином, зависящий от n переменных, - знаменатель подынтегральной формы интеграла (l);
x = [ x.,..., xn] - список переменных; Output: true, если интеграл (l) сходится абсолютно,
false - в противном случае;
1) A(P) : = NewtonPolyhedron(P, x) (используется процедура Convexhull);
2) verticesP - список вершин многогранника Ньютона A(P);
3) образовать список verticesPI, сдвинув все вершины списка verticesP на вектор I (увеличить все координаты всех точек списка verticesP на l);
4) A(Q) : = NewtonPolyhedron(Q, x) (используется программа Convexhull); {-1 (q)к,A"-1 Q)} -гиперграни A(Q);
З) {v., ... ,vd}: = NormalFan(Q, x) - веер нормалей, двойственный A(Q);
б) для каждой точки p є verticesP I и каждой
гиперграни A"- (Q) вычислить скалярное произведение вектора, отпущенного из произвольной точки гиперграни A",-1 (Q) в точку p, и вектора v^
1 < k < d ;
У) if (все скалярные произведения строго положительны) then return (true) else return(false) end if end.
Простое разбиение веера нормалей. Рассмотрим алгоритм PrimitiveNormalFan, преобразующий веер нормалей в примитивный веер нормалей в результате простого подразбиения. На входе в алгоритм имеется множество конусов X, образующих веер, на выходе - множество примитивных конусов, организующих простое подразбиение веера X.
Алгоритм состоит из двух этапов.
На первом этапе исходное множество конусов X преобразуется во множество симплициальных
конусов Е' следующим образом: пусть
C0 = (v0,..., v0m) - не симплициальный конус, т. е.
m
m > n, тогда находим вектор v = Tv1, минимизи-
i=1
руем его (сокращаем все его координаты vi на их наибольший общий делитель) и добавляем ко всем группам из (n -1) векторов множества {v0, к, v0m}, образующих конусы размерности (n - 1). Данная процедура применяется ко всем несимплициальным конусам множества X, начиная с конусов размерности 3 и заканчивая конусами размерности n. В результате мы получим множество симплициальных конусов Е' . Заметим, что так как в пространстве R2 все конусы симплициальны, то в случае интегрирования по пространству R2 первый этап автоматически опускается.
На втором этапе, пока во множестве Е' найдется конус C0 =(v0,к,v"), такой что матрица, составленная из векторов v0,..., vj , не является унимодулярной, находим точку p с целыми координатами, лежащую внутри параллелепипеда (или на его гранях, но не совпадающую ни с одной его вершиной), построенного на векторах v0, к, v" ; затем проводим из начала координат в точку p вектор v. Далее конструируем множество
Ct ^ / 1 i—1 ,+ 1
= {C,}, где C, = (v , ., v , v, v , ., v ), i = 1, ..., n; исключаем из множества C' конусы, размерность которых меньше размерности пространства; в заключение добавляем оставшиеся в C' конусы к множеству Е', а конус С0 из этого множества исключаем. В результате повторения такой процедуры получим множество примитивных конусов, которые и образуют примитивный веер.
Рассмотренный алгоритм корректен, если исходить из простых геометрических соображений.
Компьютерная реализация и примеры вычислений
Пакеты quasielliptic2 и quasielliptic3 на языке системы Maple разработаны для понижения кратности интегралов вида (1), удовлетворяющих условиям теоремы 1, по пространствам R2 и R3 соответственно. В них реализованы все алгоритмы, рассмотренные выше (за исключением алгоритма Моцкина-Бургера и алгоритма построения выпуклой оболочки, реализованногох П. А. Буров-ским [У]) в виде одноименных процедур. Все эти процедуры могут вызываться пользователем непосредственно, так как они имеют самостоятельную ценность и могут быть полезны при решении других задач.
Основной алгоритм реализован в процедуре QuasiellipticIntegration. Для ее вызова необходимо ввести следующие входные данные:
— P - полином с рациональными коэффициентами, зависящий от двух или трех переменных, - числитель коэффициента подынтегральной формы;
— Q - полином с рациональными коэффициентами, зависящий от двух или трех переменных, - знаменатель коэффициента подынтегральной формы, причем если существует полином с рациональными коэффициентами Q' и натуральное число d, такие что Q = Q'd, то вместо полинома Q в качестве исходных данных можно ввести пару [Q ', d ];
— var - список переменных, от которых зависят полиномы P и Q.
Если регулярная рациональная первообразная для подынтегральной дифференциальной формы заранее известна, то список ее коэффициентов можно ввести четвертым аргументом для ускорения работы программы.
Если исходный интеграл не удовлетворяет какому-либо из условий теоремы 1 или условию предположения 1, то выдается сообщение об ошибке и информация о том, какое из условий нарушено. Если все условия для исходного интеграла выполнены, то результатом работы процедуры для случая двух переменных является значение исходного интеграла. Для случая интегрирования по пространству R3 основной результат работы процедуры - интеграл кратности 2. Если этот интеграл вычисляется до конца стандартными средствами Maple 9, то выдается его значение; если же он до конца не вычисляется, но сводится к одномерному интегралу, то выдается именно этот интеграл. Если стандартными средствами Maple 9 двумерный интеграл не упрощается, но имеет вид (1) и удовлетворяет условиям теоремы 1, то к нему можно применить процедуру Quasiel-lipticIntegration, но уже из пакета quasielliptic2, и вычислить этот интервал до конца.
Рассмотрим примеры работы этой процедуры. Все вычисления производились в среде Maple 9 на компьютере на базе Athlon XP 1700+, 512 Mb RAM.
Сначала приведем примеры вычисления интегралов по пространству R2. Нам потребуется пакет quasielliptic2. Загрузим его следующим образом:
> read (quasielliptic2):
Пример І. Вычислим интеграл
j x2 + 4x2у2 + 4x4у4 + x4уб + xбу 4 + у2
J Л 2 2 4 4 4 4 2 21
r2 (l + х + х y + xy + xy + y)
dxdy.
Введем входные данные:
> P1 : = xЛ2+4*xЛ2*yЛ2+4*xЛ4*yЛ4+ xЛ4*yЛ6+ +xЛ6*yЛ4+yЛ2;
пі 2i/|22|/|44| 4 б і 64i 2
Pi : = x + 4xy + 4xy + x y + x y + y
> Q1 : = (1+xЛ2+xЛ2*yЛ4+xЛ4*yЛ4+ xЛ4*yЛ2+ +уЛ2)л2;
ґлл /її 2 і 2 4 і 4 4 і 42і 2\2
Q1 : = (1 + х + ху + ху + ху + У )
> уаг : = [х, у];
уаг : = [х, у]
Вызовем процедуру QuasiellipticIntegration:
> t : = tiшe():
QuasiellipticIntegration (Р1, Q1, уаг); TiшeOfCalculation : = йшє() - t;
2п
ТітеО^СаїсиїаНоп: = 0.100
Пример 2. Вычисляем интеграл
г х2 - у2 - 2х4у4+2х6у6
I т-------------------------ёхау .
я2 ( + х2 + у2 + х4у2 + х4у6 )2
Введем входные данные:
> Р2 : = хЛ2-уЛ2-2*хЛ4*уЛ4+2*хЛ6*уЛ6;
2 2 4 4 і /■> 6 6
Р2 : = х - у - 2ху + 2ху
> Q2 : = (1+хЛ2+уЛ2+хЛ4*уЛ2+хЛ4*уЛ6)Л2;
/її 2 і 2 і 4 2 і 4 6\2
Q2 : = (1 + х + у + ху + ху ) Вызовем процедуру QuasiellipticIntegration:
> t : = йшє():
QuasiellipticIntegration (Р2, Q2, уаг); TiшeOfCalculation : = tiшe() - ^
П -1 Пл/2 2
ТітеО/Саїсиїаіїоп: = 0.091 Пример 3. Вычислим интеграл
ух
J
-dxdy.
r2 (1 + x3 + 2x2y2 + 3xy4 + y6)
Введем входные данные:
> P3 : =y*x;
P3 : = yx
> Q3 : = (1+xA3+2*xA2*yA2+3*x*yA4+yA6)A3;
Q3: = (1 + x3 + 2x2y2 + 3xy4 +y6)3
Вызовем процедуру QuasiellipticIntegration:
> t : = time():
QuasiellipticIntegration (P3, Q3, var);
TimeOfCalculation : = time() - t;
Error, (in QuasiellipticIntegration) differential form denominator is not quasielliptic polynomial
TimeOfCalculation: = 0.080
Процедура возвращает сообщение о том, что знаменатель подынтегральной дифференциальной формы не является квазиэллиптическим полиномом. Это действительно так, поскольку в точке (-1, -1) полином Q3 обращается в нуль.
Рассмотрим примеры работы процедуры Qua-siellipticIntegration для случая трех переменных. Загрузим пакет quasielliptic3:
> read (quasielliptic3):
Пример 4. Вычислим интеграл
4 x 2 + 3x 2 x 2 + 3x 2 x3
f ЧХ3 +JX1 Х3 +JX2X3_________________
f { ^2 ^•A'l ^-^2^•А'З •
R
2
42
42
+ X2 + X3 + Xi X3 + X2 X
1^3
Введем входные данные:
> Pi : = 4*x3Л2+3*x1Л2*x3Л2+3*x2Л2*x3Л2;
Pi: = 4x32 + 3xi2x32 +3 x22x32
> Q1 : = (1+xlл2+x2л2+x3лз+xlлз*x3л2+ +x2л3*x3л2)л2
Q1: = (1 + xl2 + x22 + x34 + xi4x32 + x24x32)2
> var : = [xl, x2, x3];
var : = [xl, x2, x3] Вызовем процедуру QuasiellipticIntegration:
> t := time():
QuasiellipticIntegration (Pl, Ql, var); TimeOfCalculation : = time() - t;
------1----- dyldyl =1 nB\ -1,—
1 + y 2 + yl4 2 14 4
TimeOfCalculation : = 11.23б
-J J-
Пример З. Вычислим интеграл
l + x43
(l + X1
2 2 4 2 2 2 2
2 + x22 + x34 + xl2 x32 + x22 x32
-dx1dx2dx3 .
\'2 і 3 і 3 2 3 /
Введем входные данные:
> P2 : = 1+x3л3;
P2: = 1 + x34
> Q2 : = (1+xlл2+x2л2+x3лз+xlл2*x3л2+
+x2Л2*x3Л2)Л2;
Q2: = (1 + xl2 + x22 + x34 + xi2x32 +x22x32)2 Вызовем процедуру QuasiellipticIntegration:
> t : = time():
QuasiellipticIntegration (P2, Q2, var); TimeOfCalculation : = time() - t;
J J
-да -да дада
= n2-J J
1 + y22 + yl + y22 yl2
dy 2dy1 =
1 + y22 + yl2 + y22 yl2
dy Idyl = n
TimeOfCalculation : = 5.338
Таким образом алгоритмы, представленные в данной статье, позволяют понижать кратность интегралов рациональных функций со знаменателем специального вида, а в некоторых случаях точно вычислять такие интегралы. Это существенным образом расширяет класс вычисленных интегралов, что может быть полезно при решении различных математических и физических задач. Все представленные выше алгоритмы реализова-
ны автором в системе компьютерной алгебры Maple 9 для двумерного и трехмерного случая. Приведенные примеры вычисленных интегралов могут использоваться в качестве тестовых для сравнения различных методов приближенного интегрирования.
Библиографический список
1. Фам, Ф. Введение в топологические исследования особенностей Ландау / Ф. Фам. М. : Мир, 1967. 184 с.
2. Хуа, Р. Гомологии и феймановские интегралы / Р. Хуа, В. Теплиц. М. : Мир, 1969. 223 с.
3. Сергеев, А. Г. Теория твисторов и классические калибровочные поля / А. Г. Сергеев. // Мо-нополи: Топологические и вариационные методы : сб. ст. (1983-1989). М. : Мир, 1989. С. 492—555.
4. Цих, А. К. Интегралы рациональных функций по пространству R” / А. К. Цих // Докл. АН СССР. 1989. Т. 307. № 6. С. 1325-1329.
5. Ермолаева, Т. О. Интегрирование рациональных функций по R” с помощью торических компактификаций и многомерных вычетов / Т. О. Ермолаева, А. К. Цих // Мат. сб. 1996. № 9. С. 45-64.
6. Черников, С. Н. Линейные неравенства / С. Н. Черников. М. : Наука, 1968. 488 с.
7. Буровский, П. А. Алгоритм Моцкина— Бургера и вычисление выпуклых оболочек точек n -мерного пространства / П. А. Буровский // Вестн. краснояр. гос. ун-та. 2005. № 1. С. 85—92. (Серия «Физико-математические науки»).
8. Jistrand, M. Cylindrical Algebraic Decomposition - an Introduction : technical rep. / M. Jistrand ; Lincoping Univ. Lincoping, Sweden, 1995. 38 р.
9. Mishra, B. Algorithmic Algebra. Texts and Monographs in Computer Science / B. Mishra. Berlin : Springer-Verlag, 1991. 416 с.
10. Кочеткова, Т. О. Вычисление интегралов некоторых рациональных функций с помощью торических компактификаций : дис. ... канд. физ.-мат. наук / Т. О. Кочеткова. Красноярск, 1998. 74 с.
J
3
R
/ СО
l
l
M. V. Burachenko
ALGORITHM OF INTEGRAL MULTIPLICITY REDUCING FOR RATIONAL ALGEBRAICALLY EXACT DIFFERENTIAL FORM WITH QUASIELLIPTIC DENOMINATOR
It is considered a method of integral multiplicity reducing for rational algebraically exact differential from with quasielliptic denominator. It is described the algorithm and presented its implementation in the algebraic computer system Maple 9.
