Алгоритм понижения кратности интеграла рациональной алгебраически точной дифференциальной формы с квазиэллиптическим знаменателем

Каждому многограннику Ньютона A можно поставить в соответствие двойственный конический полиэдр (веер) Х_А. Дадим его определение.

С каждым ковектором v е Rⁿ * связана грань А ^v многогранника А, на которой скалярное произведение < v , и > минимально на А. Два ковектора будем считать эквивалентными, если связанные с ними грани совпадают. Замыкание классов эквивалентности ковекторов, связанных с гранью А i , образует конус в R^{n *}, порожденный элементами решетки Z . Этот конус называется двойственным к грани А i. Совокупность конусов, двойственных всем граням многогранника А, образует веер Х_А, который называется двойственным к многограннику А. Мы будем рассматривать только полные вееры, для которых объединение входящих в них конусов есть все пространство R^{n *}.

Очевидно, что образующими ребер веера А_Л являются внутренние нормали vj гиперграней

А ⁽ ” 1) е A , а набор v рождает конус a v⁷ '¹_:

, ..

, ..

., vj таких нормалей по-

., vj ) в Х_А тогда и только

тогда,    когда соответствующие грани

A ⁽ ” ^{- 1)},..., A ⁽ П - ¹)   имеют непустое пересечение

to = Zl Z iriiA ^} dx _л ... a d _x .

Q ( X 1 ,..., X n )    1 ^П

A

Размерность конуса σ будет дополни-

тельной к размерности грани A

т. е.

dima = п - dimA

Конус, число ребер которого равно его размерности, назовем симплициальным. Симплици-

альный k -мерный конус a ( v 1 ,

..

примитивным,

k

П = ^ Е X _Vv l ,0 <  X v <  1

1 1 = 1

если

., vk ) называется параллелепипед

не содержит целых точек

(точек, все координаты которых целые числа), отличных от точки 0, или, что то же самое, векторы v ¹, .., v^k можно дополнить до базиса решетки Z^{n *}. При k = n это условие эквивалентно тому, что матрица, образованная векторами v ¹, .., v^k – уни-модулярна, т. е. ее определитель равен ± 1. Веер будем называть примитивным, если любой его конус симплициален и примитивен.

Пусть d – число гиперграней многогранника A( Q ) и Z_A - веер, двойственный A ( Q ) . Обозначим через v ¹ ,..., v ^d минимально целочисленные образующие ребер веера Σ_Δ. Если веер Σ_Δ не является примитивным, то можно построить примитивный веер ~ _A , состоящий из примитивных конусов, причем каждый конус ~ _A лежит внутри некоторого конуса из S _A , и если конус из Z_A - примитивный, то он не измельчается. Пусть в результате такого подразбиения к векторам v ¹, …, v^d добавились векторы v^{d+ 1}, …, v^d+r .

Для каждого вектора v^k , 1 ≤ k ≤ d , определим

матрицу Ck = ^Ck 1 ... k_n - 1 k = ( vk ¹ ,..

., vk^{n - 1}

цами которой служат векторы v^{k 1} ,

..

¹ , v^k ., v^k

) , столб-

'^{n - 1} , поро-

ждающие вместе с v^k конус в S _A . Поскольку ~ _A - примитивный веер, то все матрицы C_k унимо-дулярны, т. е. det Ck = ±1, 1 ≤ k ≤ d .

С каждой матрицей Ck свяжем биномиальную замену координат x = yC k , что в терминах коорди-

k нат (v1 ,

..

k k

., v_n ) векторов v означает следующее:

vk1      vkn 1   vk xi = У11 ... Уп 1-1  Уп ,

vk 1       ,,kn-1    ,.k vvv xn = У1 n ...Уп-1  Упп .

.

На подынтегральное выражение в (1) будем смотреть как на дифференциальную форму, и записывать ее в виде

Определение 2 . Дифференциальная форма ω называется алгебраически точной, если существует рациональная ( n – 1) – форма φ, дифференциал d φ которой равен ω . Такую форму φ будем называть рациональной первообразной для формы ω.

Определение 3 . Пусть φ – рациональная форма ^{степени} p , т. e. ф =     £ ф ( x ) dX 1 A ... A dx , ,

1< 1 <...

функции переменной x e Rn. Если знаменатели всех функций ф являются степенями квазиэл-1... p липтического полинома Q, то мы говорим, что форма φ не растет на бесконечности, когда для каждой гиперграни а" 1 многогранника A(Q) форма ф(yCi)

не имеет полюса на гиперплоскости y_n = 0.

Введем     обозначения:      |v^k = v k +... + v k ,

[/                \       „ |vk ¹1- 1            |vk ^{n -1}1- 1 > ]

G + =H У 1 ,... У п - 1 ^)e R : y 1     ..... у П - 1      <  ° r .

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

Теорема 1 . Пусть Q – квазиэллиптический полином, не обращающийся в нуль на Rⁿ , и интеграл (1) абсолютно сходится. Если дифференциальная форма (2) имеет рациональную первообразную φ, которая не растет на бесконечности, то для интеграла (1) справедлива формула

I = ( - 1 ) ⁿ d i r ( 1 + ( - 1 ) ^{vk l}) det C_k f j - J ) [ ф ( y^Ck ) ] | y n = 0 , ^{(3) k} = ¹                       ( G +     G - )

где суммирование ведется по всем ( n – 1)-мерным граням a ^n-x ( Q ) .

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

Для случая n = 2 формулу (3) можно перепи-

сать в виде

I = ² z   J M y^{Ck )]} | y ₂ = 0

k : v k - четн R

Алгоритм понижения кратности интеграла рациональной алгебраически точной дифференциальной формы с квазиэлептическим знаменателем

На основании вышеизложенного представим алгоритм понижения кратности интеграла (1). Все вспомогательные алгоритмы, выделенные латинскими буквами, будут описаны далее, после основного алгоритма.

QuasiellipticIntegration( P, Q, x )

Input: P = P ( x 1 , … xn ) – полином, зависящий от n переменных, – числитель подынтегральной формы интеграла (1);

Q = Q ( x ₁, … x_n ) – полином, зависящий от n переменных, – знаменатель подынтегральной формы интеграла (1);

• x = [ x 1 ,… xn ] – список переменных;

Output: интеграл вида (3), кратность которого ( n – 1);

• 1)    проверка того, что Q не обращается в n нуль в R :

t 1 : = PolyNot0( Q, x );

if t 1 = false then error `знаменатель обращается в нуль`;

• 2)    проверка квазиэллиптичности полинома Q :

t 2 : = Quasielliptic( Q , x );

if t 2 = false then error `знаменатель не квазиэллиптический`;

• 3)    проверка условия предположения 1:

Δ : = NewtonPolyhedron( Q , x );

• if найдется 1 ≤ I ≤ n , такой что не суще

ствует 5 е N , такого что [ о,..., 5 , ...,0 |еА then error V i )

`условие предположения 1 не выполняется`;

• 4)    проверка того, что дифференциальная форма с коэффициентом _Q алгебраически точная, и вычисление ее рациональной первообразной φ, не растущей на бесконечности:

φ : = RationalPrimitive ( P, Q , x ) φ: = = RationalPrimitive ( Q , P , x ) ;

if ф = [0,..., 0] then error 'для исходной 123

n подынтегральной формы не существует рациональной первообразной, не растущей на беско-нечности`;

• 5)    проверка абсолютной сходимости интеграла (1):

t 3 = AbsoluteConvergence( P , Q , x );

if t 3 = false then error `интеграл не сходится абсолютно`;

• 6)    Σ : = NormalFan( Q , x ) – веер в R^{n *} , двойственный многограннику Δ, конуса которого порождены минимальными целочисленными ко-векторами v^k , k = 1, … , d внутренних нормалей к ( n - 1)-мерным граням A ^{n - 1} многогранника А;

• 7)    ~ : = Primitive Normal Fan(A) - прими-

• тивный веер, полученный в результате простого подразбиения веера Σ (Q) добавлением к его образующим ковекторов vd+1..vd+r;

• 8)    каждому ковектору v^k , k = 1, …, d , сопоставим унимодулярную матрицу

C_k = ( v k 1 ,..., v^k n ^{- 1}, v^k ) , так чтобы ковекторы v^k 1 ,..., v^k n ^{- 1}, v^k образовывали конус в ~ ( Q ) ;

• 9)    с каждой матрицей Ck , 1 ≤ k ≤ d + r свяжем биномиальную замену переменных x = yC ^k ;

• 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 соответственно)

мы будем пользоваться какой-либо одной из этих процедур в зависимости от полноты требуемых данных.

Если исходный интеграл берется по пространству R ², то подключение пакета convexhull не требуется. Для построения многогранника Ньютона (это будет выпуклый многоугольник) используется стандартная Maple-процедура convexhull программного модуля simplex. Для построения веера нормалей векторы, получающиеся при обходе многоугольника против часовой стрелки, поворачиваются на угол 90⁰ (если ( x , y ) – координаты вектора-стороны многоугольника, то (– y , x ) – координаты внутренней нормали к этой стороне) и минимизируются т. е. для каждого вектора вычисляется наименьший общий делитель его координат и все координаты на него делятся.

Проверка того, что знаменатель не обращается в нуль, и квазиэллиптичность знаменателя. На знаменатель Q ( x ) подынтегральной формы интеграла (1) накладываются два условия:

Условие 1 . Полином Q ( x ) не обращается в нуль в Rⁿ .

Условие 2 . Полином Q ( x ) является квазиэллип-тическим полиномом, т. е. для любого ненулевого ковектора v е R n * срезка Qv не обращается в нуль в торе ( R \ {0}) ⁿ .

Для проверки этих условий используется так называемый алгоритм цилиндрического алгебраического разбиения, или сокращенно CAD-алгоритм (от англ. Cylindrical Algebraic Decomposition ). Этот алгоритм является инструментом для решения систем полиномиальных уравнений и неравенств в Rⁿ и заключается в построении такого разбиения пространства Rⁿ , чтобы в каждой компоненте разбиения сохранялись знаки всех полиномов исходной системы. Под разбиением пространства будем понимать конечную совокупность непересекающихся областей пространства, объединение которых равно этому пространству

Впервые такой метод был предложен в 1948 г. Тарским, но он не был реализован практически. В 1973 г. Коллинз предложил конструктивный метод. В рамках данной статьи алгоритм был реализован в виде Maple-программы CAD.

Введем некоторые определения, необходимые для формулирования основных положений метода. Все определения даются согласно работам [8, 9]. Более детальное описание алгоритма можно найти в [9].

Входными данными для алгоритма является множество полиномов Фс Q [ x ₁,..., x n ] , где Q – поле алгебраических чисел , и список пере менных [ x 1 , ..., x n ] .

CAD- алгоритм состоит из трех основных эта пов : проекции , построения базы и расширения .

Этап проекции состоит из (n – 1)- го шага . На каждом шаге строится новое множество полино мов , зависящих от числа переменных , на единицу меньшего , чем полиномы множества , полученно го на предыдущем шаге . Множества нулей поли номов , полученных на каждом шаге этапа , – это проекции особых точек множеств полиномов предыдущего шага , т . е . точек самопересечений , изолированных точек и точек , в которых каса тельные вертикальны . Переход от множества по линомов к множеству их проекций осуществляет ся при помощи оператора proj . Чтобы его запи сать , нам понадобятся некоторые условные обо значения .

Пусть F = {f, .„, /}, где f( ( x 1 , .„, x_n ) e Q [ x _b..., x n ] , 1 ≤ i ≤ r . Эти полиномы можно рассматривать как полиномы, зависящие от одной переменной xn , коэффициентами которых являются полиномы от ( n – 1)-й переменной x 1 , …, xn– 1 , т. е. fi ( x 1 , …, xn ), 1 <  i <  r , представляется в виде f ( x 1 , ..., x_n - 1 ,x_n ) = = ft,* ( ^x 1 , -, x n - 1 ) ^x d ⁺ ... ⁺ ft ,0 ( ^x 1 ,..., x n - 1 ^{)e Q [ x} 1 ,..., x n - 1 ^][ x n ^] и f_ik ( x ₁,..., x_{n - 1} ) - коэффициент i -го полинома при k -й степени переменной, xn – полином, зависящий от переменных x 1 , …, xn –1 .

k

Обозначим    fk(x 1,..., xn) = I jx„..., xn-1)xn , j=0

где 0 ≤ k ≤ d_i .

Обозначим символом D_xn оператор дифференцирования по переменной xn .

Матрицу mat(g 1, ., gk) = [g —] размера k x k, где l = 1 + max (deg (gi)), назовем матрицей, ассо-1< i < k циированной с полиномами g1, …, gk.

Для полинома g степени s по переменной xn и полинома h степени m по переменной x_n определим psc xⁿ ( g , h ) = det ( M l ) , где матрица Ml состоит из l первых столбцов матрицы mat ( x^{s -} l ^{- 1} g ,..., g , x m ^- l ^{- 1} h ,..., h ) , ассоциированной с полиномами x s ^- l ^{- 1} g ,..., g , x m ^- l ^{- 1} h ,..., h .

Тогда proj(F) = proj1(F) и proj2(F) и рго_|'з(F), где proj1 = {fk(x1,..., xn-1)11 < i < r, 0 < k < di};

pr^oj 2 = {psc x n ^{( f}k ⁽ x_x ^,..., xn ) , Dx n ^{( f}k ^{( x} 1 ^,..., xn ) )) |1 <  i <  r , 0 <  l <  k <  d i };

p^roj 3 = ^{psc m n ^{( fk} - ⁽ x 1 ^,...^, xn ) , fk ^{j (} x ! ,..., xn ) ) |1 <  i <  j <  r , 0 <  k , k <  d., 0 <  m <  min ( k , k ) < d } .

ij i          ij i

Определенный таким образом оператор proj рекурсивно применяется (n – 1) раз к исход- ному множеству полиномов, т. е. proj0(Φ) = Φ, proj1(Φ) = proj(proj0(Φ)), …, projj(Φ) = proj(projj–1(Φ)), …, projn-1(Φ) = proj(projn-2(Φ)). Последнее множество – множество полиномов, зависящих от одной переменной x1, – конечный результат этапа проекции.

На этапе построение базы строится множество пробных точек разбиения пространства R ¹. В качестве входных данных выступает множество полиномов proj ^{n –1}( F ) – это множество полиномов, зависящих от одной переменной, получившееся в результате этапа проекции, и сама переменная, от которой зависят данные полиномы. Вычисляются действительные корни полиномов множества proj ^{n –1}( F ). Эти корни в общем случае являются алгебраическими числами и могут быть представлены в виде пары полинома и интервала, на котором этот полином имеет единственный корень – данное алгебраическое число, причем концами интервала являются рациональные числа, а полином выбирается таким образом, чтобы его степень была наименьшей. Пусть у нас получилось следующее множество алгебраических чисел § 1 e ( u i , V ] ], к , § s e ( U s , V s ] , где все U j , V j e Q . Пробные точки выбираем следующим образом: α₁ = u ₁, α 2 = ξ 1 ,…, α 2 s +1 = vs + 1. Последнее множество и является результатом этапа построения базы.

Цель этапа расширения заключается в построении разбиения пространства Rⁿ . В качестве входных данных выступает множество пробных точек, полученное на этапе построения базы: base 0 = {α 1 , …, α 2 s +1 }. Этап расширения содержит ( n – 1) шагов.

На первом шаге точки множества base 0 подставляются в полиномы множества proj ^{n –2}( F ) – проекции, полученной на ( n – 2)-м шаге этапа проекции. В результате такой подстановки для каждой точки α j , 1 ≤ j ≤ 2 s + 1 получается множество полиномов, зависящих от одной переменной x 2 . К ним снова применяется этап построения базы. Пусть в результате для каждой точки α j , 1 ≤ j ≤ 2 s + 1 получилось множество пробных точек {з/,..., e j _{+ 1} } - Тогда множество пробных точек разбиения R ² организуется следующим образом: J ^{( a}i , ^ei )• K , ^{( a}i , ^eJ11 + 1 )• ^{( a}j , ^ei )- K ,               /I

[ ( a / , e Jij + i ), к , ( a js + i , e J^{s + 1} ), к , ( a js + i , в JSj + _e + , + 1 ) j ^.

На втором шаге эти точки подставляются в полиномы из множества proj ^{n –3}( F ), и аналогичная процедура повторяется. Таким образом, после выполнения ( n – 1)-го шага мы получим множество пробных точек пространства Rⁿ . Это и есть результат CAD-алгоритма.

Корректность рассмотренного CAD-алгоритма в рамках данной статьи не доказывается. Это доказательство можно найти, например, в [9].

Вернемся теперь к проверке условий 1 и 2.

Достаточным условием того, что полином не обращается в нуль в Rⁿ , является то, что все входящие в его состав мономы c _α x ₁ ^α1 K x_n ^{α n} имеют положительные коэффициенты c α и четные степени αi по каждой переменной xi , а свободный член полинома не равен нулю. Если это условие не выполняется, то для проверки условия 1 используется CAD-алгоритм. С его помощью строится множество пробных точек для этого полинома; затем вычисляются значения исходного полинома в этих точках, и если ни в одной из этих точек полином не обращается в нуль, то можно заключить, что полином не обращается в нуль в Rⁿ . Если же хотя бы в одной пробной точке значение полинома равно нулю, то условие 1 не выполнено.

Таким образом, алгоритм PolyNot0 выглядит следующим образом:

PolyNot0( Q , x )

Input: Q = Q ( x 1 , …, xn ) – полином, зависящий от n переменных;

x = [ x ₁, …, x_n ] – список переменных;

Output: true, если полином Q не обращается в нуль в Rⁿ , false – в противном случае;

if (все мономы полинома Q имеют положительные коэффициенты и четные степени по каждой переменной) and (свободный член полинома Q равен нулю) then return (true) else

• 1)    SamplePoints : = CAD({ Q }, x ) – множество пробных точек разбиения пространства Rⁿ ;

• 2)    ValueQ – множество значений полинома Q в точках множества SamplePoints;

• 3)    if 0 не содержится в множестве ValueQ then return (true) else return (false) end if

end if end.

Достаточным условием квазиэллиптичности полинома является то, что все входящие в его состав мономы c _α x ₁ ^α1 K x_n ^{α n} имеют положительные коэффициенты c _α и четные степени α _i по каждой переменной x_i . Если это условие не выполняется, то при помощи CAD-алгоритма нужно построить множество пробных точек для всех срезок этого полинома (не имеет значения, рассматривать эти полиномы вместе или по отдельности), исключить из этого множества точки, у которых хотя бы одно координата равна нулю (это будет соответствовать исключению координатных осей из Rⁿ , так как в данном случае нас интересует только тор ( R \ {0}) ⁿ ), подставить оставшиеся точки в срезки исходного полинома и убедиться, что они не обращаются в нуль ни в одной точке из оставшихся точек. Алгоритм Quasielliptic выглядит таким образом:

Quasielliptic( Q , x )

Input: Q = Q ( x 1 , …, xn ) – полином, зависящий от n переменных;

• x = [ x 1 , …, xn ] – список переменных;

Output: true, если полином Q квазиэллиптиче-ский, false – в противном случае;

if (все мономы полинома Q имеют положительные коэффициенты и четные степени по каждой переменной) then return (true) else

• 1)    Δ( Q ):=NewtonPolyhedron( Q,x ) – многогранник Ньютона полинома Q (здесь нас будут интересовать все грани многогранника Ньютона, поэтому мы будем использовать процедуру FullConvexhull);

• 2)    F = { f 1 , …, fk } – множество срезок, соответствующих граням многогранника Δ( Q );

• 3)    SamplePoints : = CAD({ f 1 , … , fk }, x ) – множество пробных точек разбиения пространства Rⁿ ;

• 4)    SamplePointsInTor – подмножество точек из SamplePoints, у которых ни одна координата не равна нулю;

• 5)    ValueQ – множество значений полиномов множества F в точках множества SamplePoints;

• 6)    if 0 не содержится в множестве ValueQ then return (true) else return (false) end if

end if end.

Вычисление рациональной первообразной, не растущей на бесконечности. Алгоритм нахождения рациональной первообразной (Rational-Primitive), не растущей на бесконечности, основывается на методе неопределенных коэффициентов Остроградского и осуществляется в три этапа:

– определяются степени мономов, входящих в состав коэффициентов первообразной;

– выписываются полиномы с неопределенными коэффициентами и найденными степенями;

– решается уравнение на коэффициенты.

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

Рассмотрим теперь все этапы алгоритма более подробно.

Допустим, что дифференциальная форма ю = P ⁽ x ^1, K , x n ) dx,dx алгебраически точна и

Q ( X 1 , . , X n ) 1

φ – ее рациональная первообразная, не растущая на бесконечности. Представим ω в виде ~ = P ^{( x} ) dx следующим образом: если найдется

Q2 k (x)

полином Q ( x ) и натуральное число к , такие что

Q ² k ( x ) = Q ( x ) , то тогда P ( x ) = P ( x ) ; если найдется полином Q ( x ) и натуральное число к >  1, такие что Q ² k ( x ) = Q ( x ) , то тогда P ( x ) = P ( x ) • Q ( x ) ; если таких Q ( x ) и k не существует, то тогда ^{Q ( x} ) = ^{Q ( x} ) , ^k = 1, P ( x ) = P ( x ) • Q ( x ) .

Рациональную первообразную φ будем искать в виде

jxtx) dxi л .[z ]... л dxn, B (x) i=1

где x = ( x ₁, …, x_n ), A_i – некоторые полиномы; B = Q^{2k 1} - квазиэллиптический полином, не обращающийся в нуль на Rⁿ .

Пусть Δ = Δ( B ) – многогранник Ньютона. Очевидно, что многогранники Δ( Q ) и Δ( B ) гомотетичны. Таким образом, если Δ( Q ) удовлетворяет условию предположения 1, т. е. содержит по ребру на каждой координатной оси в Rⁿ , то и Δ( B ) удовлетворяет этому условию. Веера нормалей, двойственные к Δ( Q ) и Δ( B ), совпадают. Пронумеруем векторы веера нормалей, двойственного А( B ) таким образом, чтобы v^z =[ 0, . ,1, . ,0 |е Rn * ,

V i J i = 1, …, n.

Обозначим через Δ⁰ множество внутренних точек многогранника А, а через A k - множество точек, лежащих в относительной внутренности ( n -1)-мерной грани А _к = А ^ ^{- 1 )}( B ), двойственной вектору v^k .

В работе [10] доказывается теорема о том, что в предположении 1 для того чтобы дифференциальная форма φ не росла на бесконечности, ее коэффициенты должны иметь вид

Ai(x) = (- 1)*xi 'LvkAk(x), I = 1, •”, n, к=n +1

где A_i ⁰ и A^k – полиномы, многогранники Ньютона которых удовлетворяют условиям а ( A⁰ ) + I_z с А ⁰ U A ^z ₀ , a ( A k ) + 1 с А ^о _к , здесь I = ( 1, . ,1 ) e Rn , I , =( 1, . ,0, . 1 )e Rn , а v k , . , v k ^-

• V i )

координаты  вектора,  двойственного грани

AT’ (B).

Исходя из этой теоремы, определим возможные степени полиномов Ai(x), I = 1, …, n, коэффициенты которых не определены. Затем эти полиномы подставляем в уравнение jn: (- 1)z-1 f B(x)iAxl - (2к -1)A (x) 5B(x)) - P(x) = 0

i=1        V       dxz                      dxz J и решим уравнение на коэффициенты.

Если в результате решения этого уравнения все неопределенные коэффициенты получились равными нулю, то рациональной первообразной φ, не растущей на бесконечности, не существует, и мы возвращаем список из n нулей; если же не все неопределенные коэффициенты равны нулю, то возвращается список коэффициентов первообразной.

Абсолютная сходимость интеграла. В [2] доказывается следующая теорема:

Теорема 2 . Если Q – квазиэллиптический полином, не обращающийся в нуль на Rⁿ , то интеграл (1) абсолютно сходится тогда и только тогда, когда I + А ( P ) сА ⁰ ( Q ) , т. е. когда сдвиг многогранника А( P ) вдоль вектора I = ( 1, „,,1 ) е R лежит во внутренности Δ⁰( Q ) многогранника Δ( Q ).

На ее основе составлен следующий алгоритм:

AbsoluteConvergence( Q , P , x )

Input: P = P ( x 1 , …, xn ) – полином, зависящий от n переменных, числитель подынтегральной формы интеграла (1);

Q = Q ( x 1 ,…,xn ) – полином, зависящий от n переменных, – знаменатель подынтегральной формы интеграла (1);

x = [ x ₁, … , x_n ] – список переменных;

Output: true, если интеграл (1) сходится абсолютно, false – в противном случае;

• 1)    Δ( P ) : = NewtonPolyhedron( P , x ) (используется процедура Convexhull);

• 2)    verticesP – список вершин многогранника Ньютона Δ( P );

• 3)    образовать список verticesP_I , сдвинув все вершины списка verticesP на вектор I (увеличить все координаты всех точек списка verticesP на 1);

• 4)    Δ( Q ) : = NewtonPolyhedron( Q , x ) (используется программа Convexhull); { а П ¹ ( Q ), к , А d ¹ ( Q ) } -гиперграни Δ( Q );

• 5)    { v 1 , … , vd }: = NormalFan( Q , x ) – веер нормалей, двойственный Δ( Q );

• 6)    для каждой точки p е verticesP_I и каждой гиперграни А ^{п -} 1 ( Q ) вычислить скалярное произведение вектора, отпущенного из произвольной точки гиперграни а ^п - ¹ ( Q ) в точку p , и вектора v^k , 1 <  k <  d ;

• 7)    if (все скалярные произведения строго положительны) then return (true) else return(false) end if

end.

Простое разбиение веера нормалей. Рассмотрим алгоритм PrimitiveNormalFan, преобразующий веер нормалей в примитивный веер нормалей в результате простого подразбиения. На входе в алгоритм имеется множество конусов Σ, образующих веер, на выходе – множество примитивных конусов, организующих простое подразбиение веера Σ.

Алгоритм состоит из двух этапов.

На первом этапе исходное множество конусов Σ преобразуется во множество симплициальных конусов

C 0 = ^{( v} 0 ,

S' следующим образом: пусть ., v m ) - не симплициальный конус, т. е.

m m > п, тогда находим вектор v = ^ vi, минимизи-i=1

руем его (сокращаем все его координаты vⁱ на их наибольший общий делитель) и добавляем ко всем группам из ( n –1) векторов множества ^{{ v} 0 , к , ^{Vm }} , образующих конусы размерности ( n – 1). Данная процедура применяется ко всем несимплициальным конусам множества Σ, начиная с конусов размерности 3 и заканчивая конусами размерности n . В результате мы получим множество симплициальных конусов S'. Заметим, что так как в пространстве R ² все конусы симплициальны, то в случае интегрирования по пространству R ² первый этап автоматически опускается.

На втором этапе, пока во множестве S' найдется конус C ₀ = ( v 0 , к , v 0n ) , такой что матрица, составленная из векторов v ¹₀, K , v ₀ ⁿ , не является унимодулярной, находим точку p с целыми координатами, лежащую внутри параллелепипеда (или на его гранях, но не совпадающую ни с одной его вершиной), построенного на векторах v ¹₀ , K , v ₀ ⁿ ; затем проводим из начала координат в точку p вектор v . Далее конструируем множество конусов C ' = { Ci } , где Ci = ( v *, ^, i'' 1 , v , vⁱ⁺ 1 , ^, vⁿ ), i = 1, .„, n ; исключаем из множества C ' конусы, размерность которых меньше размерности пространства; в заключение добавляем оставшиеся в C ' конусы к множеству S', а конус C 0 из этого множества исключаем. В результате повторения такой процедуры получим множество примитивных конусов, которые и образуют примитивный веер.

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

Компьютерная реализация и примеры вычислений

Пакеты quasielliptic2 и quasielliptic3 на языке системы Maple разработаны для понижения кратности интегралов вида (1), удовлетворяющих условиям теоремы 1, по пространствам R ² и R ³ соответственно. В них реализованы все алгоритмы, рассмотренные выше (за исключением алгоритма Моцкина–Бургера и алгоритма построения выпуклой оболочки, реализованногох П. А. Буров-ским [7]) в виде одноименных процедур. Все эти процедуры могут вызываться пользователем непосредственно, так как они имеют самостоятельную ценность и могут быть полезны при решении других задач.

Основной алгоритм реализован в процедуре QuasiellipticIntegration. Для ее вызова необходимо ввести следующие входные данные:

– P – полином с рациональными коэффициентами, зависящий от двух или трех переменных, – числитель коэффициента подынтегральной формы;

– Q – полином с рациональными коэффициентами, зависящий от двух или трех переменных, – знаменатель коэффициента подынтегральной формы, причем если существует полином с рациональными коэффициентами Q ' и натуральное число d , такие что Q = Q ' ^d , то вместо полинома Q в качестве исходных данных можно ввести пару [ Q ', d ] ;

– var – список переменных, от которых зависят полиномы P и Q .

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

Если исходный интеграл не удовлетворяет какому-либо из условий теоремы 1 или условию предположения 1, то выдается сообщение об ошибке и информация о том, какое из условий нарушено. Если все условия для исходного интеграла выполнены, то результатом работы процедуры для случая двух переменных является значение исходного интеграла. Для случая интегрирования по пространству R ³ основной результат работы процедуры – интеграл кратности 2. Если этот интеграл вычисляется до конца стандартными средствами Maple 9, то выдается его значение; если же он до конца не вычисляется, но сводится к одномерному интегралу, то выдается именно этот интеграл. Если стандартными средствами Maple 9 двумерный интеграл не упрощается, но имеет вид (1) и удовлетворяет условиям теоремы 1, то к нему можно применить процедуру Quasiel-lipticIntegration, но уже из пакета quasielliptic2, и вычислить этот интервал до конца.

Рассмотрим примеры работы этой процедуры. Все вычисления производились в среде Maple 9 на компьютере на базе Athlon XP 1700+, 512 Mb RAM.

Сначала приведем примеры вычисления интегралов по пространству R ². Нам потребуется пакет quasielliptic2. Загрузим его следующим образом:

• >    read (quasielliptic2):

Пример 1. Вычислим интеграл с x2 + 4 x2 y2 + 4 x4 y 4 + x4 y6 + x6 y 4 + y2

J Л----- 2     2 4 , 4 4     4 2     2\ dxdy '

R 2 ( 1 + x + x y + x y + x y + y )

Введем входные данные:

• >    P1 : = x^2+4*x^2*y^2+4*x^4*y^4+ x^4*y^6+ +x^6*y^4+y^2;

P 1 : = x ² + 4 x ² y ² + 4 x ⁴ y ⁴ + x ⁴ y ⁶ + x ⁶ y ⁴ + y ²

• >    Q1 : = (1+x^2+x^2*y^4+x^4*y^4+ x^4*y^2+ +y^2)^2;

Q 1 : = (1 + x ²+ x ² y ⁴+ x ⁴ y ⁴ + x ⁴ y ²+ y ²)²

• >    var : = [x, y];

var : = [ x , y ]

Вызовем процедуру QuasiellipticIntegration:

• >    t : = time():

QuasiellipticIntegration (P1, Q1, var);

TimeOfCalculation : = time() – t;

2π

TimeOfCalculation : = 0.100

Пример 2. Вычисляем интеграл

J

R ²

x ² - y ² - 2 x ⁴ y ⁴ + 2 x ⁶ y ⁶

( 1 + x ² + y ² + x 4 y ² + x*y ⁶ ) 2

dxdy .

Введем входные данные:

• >    P2 : = x^2-y^2-2*x^4*y^4+2*x^6*y^6;

P 2 : = x ² – y ² – 2 x ⁴ y ⁴ + 2 x ⁶ y ⁶

• >    Q2 : = (1+x^2+y^2+x^4*y^2+x^4*y^6)^2;

Q 2 : = (1 + x ² + y ² + x ⁴ y ² + x ⁴ y ⁶)²

Вызовем процедуру QuasiellipticIntegration:

• >    t : = time():

QuasiellipticIntegration (P2, Q2, var);

TimeOfCalculation : = time() – t;

n — nv 2

TimeOfCalculation : = 0.091

Пример 3 . Вычислим интеграл

J ₇------------ yx-----------d-dxdy .

r ² ( 1 + x ³ + 2 x ² y ² + 3 xy⁴ + y ⁶ )

Введем входные данные:

• >    P3 : = y*x;

P 3 : = yx

• >    Q3 : = (1+x^3+2*x^2*y^2+3*x*y^4+y^6)^3;

Q 3: = (1 + x ³ + 2 x ² y ² + 3 xy ⁴ + y ⁶)³

Вызовем процедуру 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 + 3 x 2 x 2 + 3 x 2 x 3        , , ,

^{3      1 3      2 3}         dx dx dx .

I

R ³

> Q1 : = (1+x1^2+x2^2+x3^4+x1^4*x3^2+ +x2^4*x3^2)^2

Q 1: = (1 + x 1² + x 2² + x 3⁴ + x 1⁴ x 3² + x 2⁴ x 3²)²

> var : = [x1, x2, x3];

var : = [ x 1 , x 2 , x 3]

Вызовем процедуру QuasiellipticIntegration: > t := time():

QuasiellipticIntegration (P1, Q1, var);

TimeOfCalculation : = time() – t;

to to

-II-

-to -to

1 + У 2 2 + У 1 4

,      1 211 )

dy 2 dy1 = —кB I , I

2 ( 4 4 J

TimeOfCalculation : = 11.236

Пример 5 . Вычислим интеграл

1 + x 4

2   2   4   22   22

ii + x 1 + x 2 + x з + X 1 x з + x 2 x з

dx ₁ dx ₂ dx ₃

Введем входные данные: > P2 : = 1+x3^4;

P 2: = 1 + x 3⁴

> Q2 : = (1+x1^2+x2^2+x3^4+x1^2*x3^2+ +x2^2*x3^2)^2;

Q 2: = (1 + x 1² + x 2² + x 3⁴ + x 1² x 3² + x 2² x 3²)²

Вызовем процедуру QuasiellipticIntegration: > t : = time():

QuasiellipticIntegration (P2, Q2, var);

TimeOfCalculation : = time() – t;

to to                       1

I I----- 2----2-----2—2 dy 2 dy 1 =

_-№-№ 1 + y 2 + y1 + y 2 y1

= п

to to

-II-

-to -to

1 + y 2 2 + y 1² + y 2² y 1²

dy 2 dy 1 = п²

TimeOfCalculation : = 5.338

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

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