1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 63, № 1, 2023

Журнал вычислительной математики и математической физики, 2023, T. 63, № 1, стр. 93-101

Аналитическое исследование кубатурных формул на сфере в системах компьютерной алгебры

Р. Э. Байрамов 1, Ю. А. Блинков 1, 2, И. В. Левичев 1, М. Д. Малых 1, 3, *, В. С. Мележик 3

1 РУДН
117198 Москва, ул. Миклухо-Маклая, 6, Россия

2 Саратовский гос. ун-т им. Н.Г. Чернышевского
410012 Саратов, ул. Астраханская, 83, Россия

3 ОИЯИ
141980 М.о., Дубна, Россия

* E-mail: malykh_md@pfur.ru

Поступила в редакцию 24.04.2022
После доработки 24.04.2022
Принята к публикации 10.09.2022

Полный текст (PDF)

Аннотация

Задача об отыскании весов и узлов кубартурных формул заданного порядка на единичной сфере, инвариантных относительно групп вращения икосаэдра (задача А.С. Попова) исследуется аналитически в системах компьютерной алгебры. Алгоритм Попова сведения задачи к системе нелинейных уравнений реализован в известной системе компьютерной алгебры Sage. Показано, что в Sage трудности с исследованием полученной системы нелинейных алгебраических уравнений возникают, начиная с порядка аппроксимации, равного 23. Показано также, что задача Попова при таком порядке приводит к полиномиальному идеалу, базис Грёбнера для которого содержит многочлены с экстремально большими целыми коэффициентами, что делает ее весьма трудной для исследования стандартными инструментами, реализованными в Sage. Этот базис найден в нашей системе компьютерной алгебры – GInv, новая версия которой была передана в общественный доступ одним из авторов статьи в 2021 г. Это позволило далее полностью описать множество решений задачи Попова в Sage. Проведено сравнение найденных нами точных решений с решениями, найденными Поповым численно. Обсужден потенциал использования задачи Попова как тестовой задачи для систем, специализирующихся на вычислении базиса Грёбнера. Библ. 19.

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

1. ВВЕДЕНИЕ

Задача об отыскании кубатурных формул, инвариантных относительно групп симметрии правильных многогранников, была рассмотрена в серии статей А.С. Попова, где были описаны алгоритмы отыскания кубатурных формул и их численная реализация для наиболее интересных групп, а именно, для группы вращений октаэдра с и без инверсии (см. [1], [2]) и группы вращений икосаэдра с и без инверсии (см. [3], [4]). Эти результаты были подытожены автором в [5]. Поэтому задачу об отыскании кубатурных формул заданного порядка, инвариантных относительно групп симметрий многогранников, мы будем называть задачей Попова. Эта задача была также предметом недавних исследований М. Грефа (см. [6]–[8]).

Напомним кратко, в чем состоит проблема. Пусть $S$ – единичная сфера

${{x}^{2}} + {{y}^{2}} + {{z}^{2}} = 1.$
Рассмотрим интеграл
$U(f) = \frac{1}{{4\pi }}\iint\limits_S \,fds,$
где $ds$ – элемент поверхности единичной сферы. Под кубатурной формулой для этого интеграла понимают выражение вида
$V(f) = \sum\limits_{i = 1}^N \,{{A}_{i}}f({{s}_{i}}),$
которое можно использовать для приближенного вычисления интеграла $U(f)$. Здесь ${{s}_{i}}$ – точки на сфере $S$, именуемые узлами, а ${{A}_{i}}$ – положительные числа, именуемые весами. Говорят, что кубатурная формула $V$ имеет алгебраический порядок точности $n$, если равенство
$U(f) = V(f)$
верно для всех многочленов $f \in \mathbb{Q}[x,y,z]$, степень которых $\partial f \leqslant n$. Говорят, что кубатурная формула $V$ инвариантна относительно группы симметрий $\mathfrak{G}$ сферы, если относительно этой группы инвариантно множество ее узлов и

$\sum\limits_{i = 1}^N \,{{A}_{i}}f(T{{s}_{i}}) = \sum\limits_{i = 1}^N \,{{A}_{i}}f({{s}_{i}}),\quad T \in \mathfrak{G}.$

Хотя исторически кубартурные формулы вышли из задачи о приближенном вычислении интеграла $U(f)$ на сфере, их приложения много шире. Дело в том, что численное исследование дифференциальных уравнений в той или иной двумерной области опирается на возможность введения на ней сетки (mesh), узлы которой распределены более менее равномерно. Однако, поскольку число правильных многогранников не бесконечно, введение на сфере такой равномерной сетки представляет известные сложности. Формулы, найденные Поповым, дают такое распределение узлов и активно используются в работах В.С. Мележика (см. [9]–[11]) для построения двумерного DVR базиса, который хорошо зарекомендовал себя в аппроксимации трехмерного стационарного уравнения Шрёдингера, моделирующего спектр атома водорода в сильных скрещенных электрическом и магнитных полях, и нестационарного уравнения Шрёдингера с тремя пространственными переменными, моделирующего генерацию высоких гармоник атомом водорода в сильном эллиптически поляризованном лазерном поле.

В работах А.С. Попова решение задачи Попова сводится к исследованию системы нелинейных алгебраических уравнений, алгоритм составления которой описан в его статьях. Затем эта задача решалась численно на вычислительной технике Сибирского суперкомпьютерного центра (см. [5]). М. Греф (см. [8]) тоже занимался численным исследованием возникающих систем нелинейных уравнений. Современные системы компьютерной алгебры, вообще говоря, в состоянии предложить богатый инструментарий для составления систем алгебраических уравнений и исследования множества их решения в различных полях. В настоящей работе задача Попова исследуется в доступных и свободных системах компьютерной алгебры: популярной системе общего назначения Sage (https://www.sagemath.org) и созданной совсем недавно системе GInv (https://github.com/blinkovua/GInv), ориентированной на исследование полиномиальных идеалов.

Исследования системы нелинейных уравнений в современных системах компьютерной алгебры требуют на вычисления базисов Грёбнера идеалов, порожденных левыми частями уравнений этих систем (см. [12]). Для вычисления базиса Грёбнера в Sage используется реализация алгоритма Бухбергера (см. [12]), пришедшая в эту систему из системы Singular. Алгоритм Бухбергера – старейший, базовая версия этого алгоритма оставляет много свободы в осуществлении вычислительного процесса, значительные улучшения получаются за счет реализации критериев уменьшения числа S-многочленов, подлежащих рассмотрению. В конце 1990-х гг. в качестве альтернативы алгоритму Бухбергера были предложены инволютивные алгоритмы (см. [13]–[16]), которые реализованы в системе GInv. По инициативе В.П. Гердта и Ю.А. Блинкова в 2005 г. был основан проект GInv (http://invo.jinr.ru), в рамках которого было разработано программное обеспечение для вычисления инволютивных базисов, написанное в виде С++ модуля для языка Python. Недавно эта система была существенно переработана одним из авторов настоящей статьи, новая версия находится в открытом доступе по адресу https://github.com/blinkovua/GInv. В новую версию системы GInv добавлены механизмы динамического перераспределения памяти, позволяющие существенно, до нескольких раз, ускорить вычисления. Предлагаемый подход принципиально отличается от других алгоритмов известных под общим названием “сборка мусора” (garbage collection) (см. [17], [18]) и основан на реализации объектно-ориентированного программирования в C++.

Разумеется, мы сравнили базисы Грёбнера, полученные в Sage и GInv для задачи Попова при $n = 19$, и убедились в их совпадении с точностью до редукции лишних базисных элементов. Этот тест, наряду со стандартными “синтетическими” тестами, убедил нас в корректности работы новой системы.

В теории мы можем вычислить базис Грёбнера любой системы за конечное число шагов, однако на практике на эти шаги очень часто не хватает ресурсов, в первую очередь объема памяти компьютера. В настоящее время десяток неизвестных дает ту границу, после пересечения которой нельзя быть уверенным в том, что ресурсов современных PC хватит. Задача Попова при порядке аппроксимации, равном 20, как раз лежит на этой границе, что придает исследованию задачи Попова в системах компьютерной алгебры особый интерес.

2. КУБАТУРНЫЕ ФОРМУЛЫ, ИНВАРИАНТНЫЕ ОТНОСИТЕЛЬНО ГРУППЫ ВРАЩЕНИЙ ИКОСАЭДРА

Для отыскания кубатуры, инвариантной относительно группы вращений икосаэдра, Попов (см. [3]) зафиксировал на сфере один икосаэдр, связал с ним первые точки сети, а кубатуру переписал в виде

(1)
$V(f) = {{A}_{0}}\sum\limits_{j = 1}^{12} \,f({{a}_{{0j}}}) + {{B}_{0}}\sum\limits_{j = 1}^{20} \,f({{b}_{{0j}}}) + {{C}_{0}}\sum\limits_{j = 1}^{30} \,f({{c}_{{0j}}}) + \sum\limits_{i = 1}^M \,{{A}_{i}}{\kern 1pt} \sum\limits_{j = 1}^{60} \,f({{a}_{{ij}}}),$
где 12 точек ${{a}_{{0j}}}$ лежат в вершинах вписанного в сферу икосаэдра, 20 точек ${{b}_{{0j}}}$ отвечают центрам граней икосаэдра, 30 точек ${{c}_{{0j}}}$ отвечают серединам ребер икосаэдра, а 60 точек ${{a}_{{ij}}}$ отвечают точкам общего положения на гранях икосаэдра, точки ${{a}_{{ij}}}$ получаются из ${{a}_{{i1}}}$ путем вращения икосаэдра. При таком выборе условие инвариатности квадратурной формулы относительно группы вращений икосаэра выполняется автоматически.

В силу теоремы Соболева (см. [3, p. 435]), для того чтобы кубатурная формула (1) имела порядок $n$, необходимо и достаточно, чтобы она была точна для всех инвариантных относительно группы вращения многочленов степени не выше $n$. Опираясь на результаты Клейна (см. [19]), все эти инвариантные многочлены можно описать как многочлены относительно трех многочленов, которые в работе Попова обозначены как $u,\;{v}$ и $w$. Мы не будем выписывать здесь громоздкие выражения для этих многочленов, отметим лишь, что многочлены эти не независимы, но связаны алгебраическим уравнением

(2)
${{w}^{2}} = - 16{{v}^{2}} + 40{{u}^{2}}v - 25{{u}^{4}} + 20u{{v}^{2}} - 45{{u}^{3}}v + 27{{u}^{5}} - {{v}^{3}}.$
Поскольку $u({{a}_{{ij}}})$ не зависит от $j$, мы можем обозначить это значение как ${{u}_{i}}$. Приняв аналогичные обозначения для $v$ и $w$, заключаем: кубатурная формула имеет порядок $n$ тогда и только тогда, когда
$V({{u}^{p}}{{v}^{q}}{{w}^{r}}) = 12{{A}_{0}}u_{0}^{p}v_{0}^{q}w_{0}^{r} + \ldots + 60\sum\limits_{i = 1}^M \,{{A}_{i}}u_{i}^{p}v_{i}^{q}w_{i}^{r}$
для всех выражений ${{u}^{p}}{{v}^{q}}{{w}^{r}}$, степень которых относительно $x,\;y,\;z$ не превосходит $n$. Левые части этих равенств нетрудно вычислить явно. Поэтому получается система алгебраических уравнений относительно весов и значений инвариантных многочленов в узлах сети. Таким образом, для их отыскания получается система нелинейных алгебраических уравнений, которую мы далее будем называть системой Попова.

На основе этих соображений Попов разработал и описал в [3] алгоритм, при реализации которого в системе Sage мы не столкнулись с трудностями, достойными упоминания. В [3] приведены результаты вычислений для порядков $n$ с 19 до 25 в виде таблиц с десятичными дробями, случай $n = 19$ исследован аналитически.

2.1. Кубатура при $n = 19$

Кубатура $n = 19$, $M = 2$, ${{B}_{0}} = {{C}_{0}} = 0$ была исследована Поповым аналитически. Наша функция popoff_eqs(n) возвращает список, первым элементом которого служит список всех искомых величин (весов и значений инвариантов), затем следует список всех уравнений системы. Для $n = 19$ мы имеем:

sage: load(“popoff.sage”)

sage: popoff_eqs(19)

[[A0, A1, A2, u1, u2, v1, v2, w1, w2],

[12*A0 + 60*A1 + 60*A2 – 1,

60*A1*w1 + 60*A2*w2,

60*A1*v1 + 60*A2*v2 – 256/231,

60*A1*u1 + 60*A2*u2 – 16/21,

60*A1*u1*v1 + 60*A2*u2*v2 – 18432/17017,

60*A1*u1^2 + 60*A2*u2^2 – 2048/3003,

60*A1*u1^3 + 60*A2*u2^3 – 8192/12597,

27*u1^5 – 25*u1^4 – 45*u1^3*v1 + 40*u1^2*v1 + 20*u1*v1^2 – v1^3 – 16*v1^2 – w1^2,

27*u2^5 – 25*u2^4 – 45*u2^3*v2 + 40*u2^2*v2 + 20*u2*v2^2 – v2^3 – 16*v2^2 – w2^2]]

В Sage базис Грёбнера идеала $J$, порожденный этими уравнениями (deglex-порядок на мономах принят по умолчанию), вычисляется без каких-либо затруднений и существенных затрат времени:

sage: eqs=popoff_eqs(19)

sage: K=PolynomialRing(QQ,eqs[0])

sage: J=K*eqs [1]

sage: J.groebner_basis()

[v2^2 – 3038016/1282633*v2 + 197738496/239852371,

v2*w1 – 22131456/1115081*w1 – 20793600/1115081*w2,

w1^2 – 282712985600/93773402523273*v2 – 17460378265190400/331341213770339443,

v2*w2 + 20793600/1115081*w1 + 12153087936/695303689*w2,

w1*w2 + 910295040000/16191714009097,

w2^2 + 282712985600/93773402523273*v2 – 349759628902400/5845208757284017,

A0 – 125/19656,

A1 + 85936411/416549952000*v2 – 2153726/271191375,

A2 – 85936411/416549952000*v2 – 5389889/723177000,

u1 + 19/60*v2 – 403344/337535,

u2 – 19/60*v2 – 416/935,

v1 + v2 – 3038016/1282633]

Отметим, что в данном случае “упрощение” системы путем перехода к целым коэффициентам не только не требуется, но и несколько затягивает вычисления. Интересно отметить и огромные коэффициенты, которые появляются в базисе уже при $n = 19$.

Быстрота вычисления базиса Грёбнера позволяет без хлопот вычислить размерность идеала, которая равна нулю. После этого без ошибок можно применить к идеалу метод variety и получить множество решений системы уравнений Попова над полем вещественных алгебраических чисел $\mathbb{A}$ в виде списка L:

sage: L=J.variety(AA)

Таким образом, при $n = 19$ система Попова имеет четыре решения. У всех четырех вес ${{A}_{0}}$ имеет одно и то же рациональное значение

${{A}_{0}} = \frac{{125}}{{19656}},$
указанное Поповым. Это число нетрудно найти в Sage по минимальному многочлену для алгебраического числа ${{A}_{0}}$:

sage: L[0][A0].minpoly()

x – 125/19656

Два других веса ${{A}_{1}}$ и ${{A}_{2}}$ являются корнями одного и того же квадратного уравнения:

sage: L[0][A1].minpoly()

x^2 – 1513/98280*x + 42093683/710738355600

sage: L[0][A2].minpoly()

x^2 – 1513/98280*x + 42093683/710738355600

Поэтому во всех четырех решениях веса имеют одни и те же значения с точностью до перестановки ${{A}_{1}}$ и ${{A}_{2}}$. Эти значения указаны в таблице Попова как одно решение, сами значения совпадают до последнего знака, указанного в таблице Попова.

В данном случае кубатура (1) имеет вид

$V(f) = {{A}_{0}}\sum\limits_{j = 1}^{12} \,f({{a}_{{0j}}}) + {{A}_{1}}\sum\limits_{j = 1}^{60} \,f({{a}_{{ij}}}) + {{A}_{2}}\sum\limits_{j = 1}^{60} \,f({{a}_{{ij}}}).$
Поэтому из одного решения системы можно получить второе путем одновременной перестановки весов ${{A}_{1}}$ и ${{A}_{2}}$ и значений инвариантов ${{u}_{1}},\;{{v}_{1}},\;{{w}_{1}}$ и ${{u}_{2}},\;{{{v}}_{2}},\;{{w}_{2}}$. Куда менее очевидно, что из одного решения системы можно получить второе путем одновременного изменения знаков у ${{w}_{1}}$ и ${{w}_{2}}$. Это преобразование соответствует отражению на сфере, т.е. получающиеся при этом узлы будут совпадать с точностью до отражения. В итоге получается четыре различных решения системы Попова, которые в таблице Попова отмечены как одно.

2.2. Кубатура при $n = 20$

При $n = 20$ мы имеем $M = 2$ и ${{C}_{0}} = 0$. Система Попова имеет вид

sage: popoff_eqs(20)

[[A0, B0, A1, A2, u1, u2, v1, v2, w1, w2],

[12*A0 + 60*A1 + 60*A2 + 20*B0 – 1,

60*A1*w1 + 60*A2*w2, 60*A1*v1 + 60*A2*v2 + 5120/81*B0 – 256/231,

60*A1*v1^2 + 60*A2*v2^2 + 1310720/6561*B0 – 1900544/969969,

60*A1*u1 + 60*A2*u2 + 640/27*B0 – 16/21,

60*A1*u1*v1 + 60*A2*u2*v2 + 163840/2187*B0 – 18432/17017,

60*A1*u1^2 + 60*A2*u2^2 + 20480/729*B0 – 2048/3003,

60*A1*u1^3 + 60*A2*u2^3 + 655360/19683*B0 – 8192/12597,

27*u1^5 – 25*u1^4 – 45*u1^3*v1 + 40*u1^2*v1 + 20*u1*v1^2 – v1^3 – 16*v1^2 – w1^2,

27*u2^5 – 25*u2^4 – 45*u2^3*v2 + 40*u2^2*v2 + 20*u2*v2^2 – v2^3 – 16*v2^2 – w2^2]]

Над полем вещественных алгебраических чисел $\mathbb{A}$ эта система имеет четыре решения. У всех четырех вес ${{A}_{0}}$ имеет одно и то же значение, которое совпадает с указанным в таблице Попова, но, в отличие от случая $n = 19$, не является рациональным числом:

sage: L[0][A0].minpoly()

x^2 – 7375/453024*x + 1328125/20777492736

На самом деле, мы ожидали, что старший вес ${{A}_{0}}$ всегда является рациональным числом, поскольку в сумме (1) он занимает особое положение, в то время как младшие веса ${{A}_{1}}$ и ${{A}_{2}}$ можно менять местами.

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

Два других веса ${{A}_{1}}$ и ${{A}_{2}}$ являются корнями одного и того же уравнения 4-й степени:

sage: L[0][A1].minpoly()

x^4 – 11885261/482017536*x^3 + 202915177968491/1033960700403671040*x^2 – 2389552908217687/4766422781400344064000*x – 26951199206199859/610712218135263284232192000

sage: L[0][A2].minpoly() x^4

– 11885261/482017536*x^3 + 202915177968491/1033960700403671040*x^2 – 2389552908217687/4766422781400344064000*x – 26951199206199859/610712218135263284232192000

Эти уравнения имеют четыре вещественных корня, но лишь два из них дают вещественные значения для инвариантов. Структура решения над $\mathbb{A}$ та же, что и при $n = 19$: из одного решения получаются другие путем перестановки индексов $1$ и $2$ и путем одновременной смены знаков у ${{w}_{1}}$ и ${{w}_{2}}$.

2.3. Кубатура при $n = 23$

Обратимся теперь к случаю $n = 23$, когда $M = 3$, ${{B}_{0}} = {{C}_{0}} = 0$, как и при $n = 19$. Мы ожидали, что этот случай должен быть в чем-то схож с простейшим случаем $n = 19$:

sage: load(“popoff.sage”)

sage: popoff_eqs(23)

[[A0, A1, A2, A3, u1, u2, u3, v1, v2, v3, w1, w2, w3],

[12*A0 + 60*A1 + 60*A2 + 60*A3 – 1,

60*A1*w1 + 60*A2*w2 + 60*A3*w3,

60*A1*v1 + 60*A2*v2 + 60*A3*v3 – 256/231,

60*A1*v1^2 + 60*A2*v2^2 + 60*A3*v3^2 – 1900544/969969,

60*A1*u1 + 60*A2*u2 + 60*A3*u3 – 16/21,

60*A1*u1*w1 + 60*A2*u2*w2 + 60*A3*u3*w3,

60*A1*u1*v1 + 60*A2*u2*v2 + 60*A3*u3*v3 – 18432/17017,

60*A1*u1^2 + 60*A2*u2^2 + 60*A3*u3^2 – 2048/3003,

60*A1*u1^2*v1 + 60*A2*u2^2*v2 + 60*A3*u3^2*v3 – 8126464/7436429,

60*A1*u1^3 + 60*A2*u2^3 + 60*A3*u3^3 – 8192/12597,

27*u1^5 – 25*u1^4 – 45*u1^3*v1 + 40*u1^2*v1 + 20*u1*v1^2 – v1^3 – 16*v1^2 – w1^2,

27*u2^5 – 25*u2^4 – 45*u2^3*v2 + 40*u2^2*v2 + 20*u2*v2^2 – v2^3 – 16*v2^2 – w2^2,

27*u3^5 – 25*u3^4 – 45*u3^3*v3 + 40*u3^2*v3 + 20*u3*v3^2 – v3^3 – 16*v3^2 – w3^2]]

В системе Sage нам не удалось найти базис Грёбнера этой системы. Эффективность алгоритмов отыскания базиса существенно зависит от порядка следования неизвестных. Мы пробовали различные варианты, но, конечно, не перебрали все 13(!) вариантов.

Процессы отыскания базиса Грёбнера при deglex порядке на мономах и

${{A}_{0}} > {{A}_{1}} > {{A}_{2}} > {{A}_{3}} > {{u}_{1}} > {{u}_{2}} > {{u}_{3}} > {{v}_{1}} > {{v}_{2}} > {{v}_{3}} > {{w}_{1}} > {{w}_{2}} > {{w}_{3}}$
мы запустили параллельно в двух системах Sage и GInv на компьютере “Научного центра вычислительных методов в прикладной математике” РУДН. С задачей успешно справился только GInv, вернув нам базис $B$, состоящий из 107 многочленов небольшой степени (degree), но с огромными целыми коэффициентам порядка ${{10}^{{161}}}$. Среди них два многочлена были линейными, 42 многочлена – квадратичными и 63 многочлена – кубическими.

Поскольку на данном этапе система GInv не имеет разработанного инструментария для дальнейшего использования найденного базиса, мы передали его в Sage и определили идеал $J$. Работа с идеалом $J$, порожденным 107 многочленами с 13 переменными и огромными целыми коэффициентами, является экстремально трудной. Часто применение к этому идеалу стандартных методов приводило к долгим вычислениям, которые мы прерывали через несколько минут. Тем не менее несколько методов позволили полностью решить задачу Попова.

Прежде всего мы решили проверить, является ли найденный базис базисом Грёбнера с помощью стандартного метода basis_is_groebner, однако за несколько минут ответа не получили и прервали процесс. Затем мы решили вычислить базис Грёбнера идеала $J$ при том же мономиальном порядке, что и в GInv, с помощью метода groebner_basis и достигли успеха буквально за несколько секунд:

sage: J.groebner_basis()

Polynomial Sequence with 83 Polynomials in 13 Variables

В Sage по умолчанию выводится редуцированный базис, поэтому уменьшения числа элементов и следовало ожидать.

Применение к идеалу метода variety вновь приводит к излишне затратным вычислениям, которые были прерваны. Однако нам удалось обойти эту трудность, поскольку в Sage за несколько секунд нам удалось вычислить базис Грёбнера относительно lex порядка, т.е. фактически перейти от базиса Грёбнера относительно deglex мономиального порядка, вычисленному в GInv, к базису Грёбнера относительно lex порядка. Мы взяли такой порядок, чтобы ${{A}_{0}}$ был младшим:

sage: K=PolynomialRing(QQ, [w1, w2, w3, u1, u2, u3, v1, v2, v3, A1, A2, A3, A0], order='lex')

sage: Jlex=Klex*B

Базис Грёбнера при таком выборе мономиального порядка имеет всего 13 элементов:

sage: Blex=Jlex.groebner_basis()

sage: len(Blex)

13

Напомним, что неизвестных у нас тоже 13.

Базис Грёбнера при lex порядке представляет собой систему многочленов , приведенную к треугольному виду. Разберем этот базис с конца. Последней элемент базиса дает многочлен 3‑й степени для ${{A}_{0}}$:

sage: Blex[-1]

A0^3 – 453214375/25273614366*A0^2

+ 6016015625/57483863813376*A0

– 2090576171875/10608532099030914048

Это многочлен имеет три вещественных корня:

sage: QQ[A0](Blex[-1]).roots(AA)

[(0.004164880157580243?, 1),

(0.006620018869643640?, 1),

(0.007147414429679166?, 1)]

Один из них, наименьший, совпадает с точностью до последнего знака со значением веса ${{A}_{0}}$ в таблице Попова.

Следующий многочлен

sage: Blex[-2]

A3^3 + 1/5*A3^2*A0 – 1/60*A3^2 – …

содержит ${{A}_{3}}$ и ${{A}_{0}}$ и позволяет для любого из трех значений ${{A}_{0}}$ найти по три значения ${{A}_{3}}$. При всех трех возможных значениях ${{A}_{0}}$ этот многочлен относительно ${{A}_{3}}$ имеет вещественные корни:

sage: [(A01,a),(A02,b),(A03,c)]=QQ[A0](Blex[-1]).roots(AA)

sage: AA[A3](Blex[-2].subs(A0 1)).roots()

[(0.0050776793060011?, 1),

(0.0053495468022166?, 1),

(0.0054064645269329?, 1)]

sage: AA[A3](Blex[-2].subs(A0 2)).roots()

[(0.001590837633796970?, 1),

(0.006614776888519542?, 1),

(0.007137048370421428?, 1)]

sage: AA[A3](Blex[-2].subs(A0 3)).roots()

[(0.0002487875296393600?, 1),

(0.007210715954234012?, 1),

(0.007777680296857461?, 1)]

Вообще говоря, далее следовало бы рассматривать девять случаев и из следующих двух уравнений найти ${{A}_{2}},\;{{A}_{1}}$. Однако задача Попова имеет естественную симметрию: если взять $i,j > 0$ и на решении переставить

${{A}_{i}},{{u}_{i}},{{v}_{i}},{{w}_{i}}$
и
${{A}_{j}},{{u}_{j}},{{v}_{j}},{{w}_{j}},$
то значение суммы (1) не изменится, т.е. получится опять решение задачи Попова. Как отмечалось выше, в таблице Попова указывается одно решение по модулю этих перестановок. В данном случае эта симметрия ведет к тому, что веса ${{A}_{1}}$ и ${{A}_{2}}$ удовлетворяют тому же уравнению 3-й степени, что и ${{A}_{3}}$, причем ${{A}_{1}},\;{{A}_{2}}$ и ${{A}_{3}}$ являются тремя различными корнями уравнения
(3)
${{A}^{3}} + \frac{1}{5}{{A}_{0}}{{A}^{2}} - \frac{1}{{60}}{{A}^{2}} - \ldots = 0,$
где в качестве ${{A}_{0}}$ может быть один из трех корней уравнения

(4)
$A_{0}^{3} - \frac{{453214375}}{{25273614366}}A_{0}^{2} + \frac{{6016015625}}{{57483863813376}}{{A}_{0}} - \frac{{2090576171875}}{{10608532099030914048}} = 0.$

Следующие шесть многочленов, с пятого до десятого с конца, – линейные относительно ${{u}_{1}},\;{{u}_{2}},\;{{u}_{3}}$, ${{{v}}_{1}},\;{{{v}}_{2}},\;{{{v}}_{3}}$, они позволяют выразить эти шесть переменных через найденные выше веса.

Одиннадцатый с конца многочлен дает уравнение вида

$w_{3}^{2} = f({{A}_{0}},{{A}_{3}}).$
Имеет ли оно вещественный корень или нет, зависит от знака многочлена $f({{A}_{0}},{{A}_{3}})$. Подставляя вместо ${{A}_{0}}$ и ${{A}_{3}}$ корни выписанных выше уравнений, мы видим, что вещественное значение для ${{w}_{3}}$ получается в двух случаях: когда в качестве ${{A}_{0}}$ берется наименьший или наибольший корень. При среднем корне получается мнимое значение.

Первые два многочлена базиса – линейные относительно ${{w}_{1}}$ и ${{w}_{2}}$ соответственно, из них можно однозначно и без введения комплексных чисел найти ${{w}_{1}}$ и ${{w}_{2}}$.

Таким образом, над $\mathbb{A}$ множество решений устроено следующим образом. Вес ${{A}_{0}}$ равен наибольшему или наименьшему из трех корней уравнения 3-й степени (4). Веса ${{A}_{1}},\;{{A}_{2}},\;{{A}_{3}}$ являются тремя корнями уравнения (3). По ним однозначно определяются вещественные значения инвариантов $u$ и $v$, а значения $w$ – с точностью до знака. Всего, с точностью до перестановок, два корня: корень с

${{A}_{0}} = 0.004164880157580243,$
указанный в таблице Попова, и корень
${{A}_{0}} = 0.007147414429679166,$
в таблице отсутствующий. Таким образом, в случае $n = 23$ задача Попова допускает два решения, различающихся весом первого члена. Численно из этих решений было найдено лишь одно, отвечающее наименьшему весу.

3. ЗАКЛЮЧЕНИЕ

Задача Попова дает в наше распоряжение набор систем нелинейных алгебраических уравнений, сложность решения которых растет вместе с ростом порядка аппроксимации $n$. Мы имеем численные решения этой системы в интервале от $n = 19$ до $n = 25$. Случай $n = 19$ был исследован аналитически Поповым на руках, случай $n = 20$ успешно исследован в системе Sage. При этом, вопреки нашим ожиданиям, старший вес ${{A}_{0}}$ оказался иррациональным числом. Сравнение аналитики и численных расчетов подтверждает, что численные расчеты Попова выполнены с очень высокой точностью. Уравнения Попова при $n = 23$ являются слишком сложными для того, чтобы их можно было исследовать стандартными средствами Sage/Singular, зато с ней справляется инволюционный алгоритм, реализованный в GInv. В нашем распоряжении имеются численные решения, поэтому задача Попова является превосходным тестом, позволяющим оценить эффективность различных алгоритмов вычисления базиса Грёбнера.

Применение инволютивных алгоритмов позволило описать все множество решений задачи Попова при $n = 23$ и, в частности, найти ранее неизвестное решение.

Следует особо подчеркнуть, что этот тест не является искусственным, специально придуманным для тестирования систем компьютерной алгебры, а результаты сами по себе представляют научный интерес. Нельзя не отметить значительное отличие параметров, характеризующих вычисление базиса Грёбнера, при исследовании задачи Попова и стандартных тестов (http://invo.jinr.ru/ginv/benchmark.html): коэффициенты базисных многочленов оказываются огромными, а степени, напротив, остаются небольшими.

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

  1. Popov A.S. The search for the sphere of the best cubature formulae invariant under octahedral group of rotations // Siberian J. of Numer. Math. / Sib. Branch of Russ. Acad. of Sci. 2002. V. 5. № 4. P. 367–0372.

  2. Popov A.S. The search for the best cubature formulae invariant under the octahedral group of rotations with inversion for a sphere // Siberian J. of Numer. Math. / Sib. Branch of Russ. Acad. of Sci. 2005. V. 8. № 2. P. 143–148.

  3. Popov A.S. Cubature formulas on a sphere invariant under the icosahedral rotation group // Numer. Analys. A-ppl. 2008. V. 1. P. 355–361.

  4. Popov A.S. Cubature Formulas Invariant under the Icosahedral Group of Rotations with Inversion on a Sphere // Numer. Analys. Appl. 2017. V. 10. P. 339–346.

  5. Popov A.S. Cubature formulas on a sphere invariant under the symmetry groups of regular polyhedrons // Siberian Electron. Math. Rep. 2017. V. 14. P. 190–198.

  6. Gräf M., Potts D. Sampling sets and quadrature formulae on the rotation group // Numer. Funct. Anal. and Optimizat. 2009. V. 30. № 7–8. P. 665–688.

  7. Gräf M., Kunis St., Potts D. On the computation of nonnegative quadrature weights on the sphere // Appl. and Comput. Harmonic Anal. 2009. V. 27. № 1. P. 124–132.

  8. Gräf M. Efficient Algorithms for the computation of Optimal Quadrature Points on Riemannian Manifolds : Ph. D. thesis / Manuel Gräf; Fakultät für Mathematik der Technischen Universität Chemnitz. Chemnitz, 2013. P. 2.

  9. Melezhik V.S. New method for solving multidimensional scattering problem // J. of Comput. Phys. 1991. V. 92. № 1. P. 67–81. Access mode: https://www.sciencedirect.com/science/article/pii/002199919190292S.

  10. Shadmehri S., Saeidian Sh., Melezhik V.S. 2D nondirect product discrete variable representation for Schrödinger equation with nonseparable angular variables // J. of Phys. B: Atomic, Molecular and Optic. Phys. 2020. V. 53. № 8. P. 085001.

  11. Melezhik V.S. Improving efficiency of sympathetic cooling in atom- ion and atom-atom confined collisions // Phys. Rev. A. 2021. May. V. 103. P. 053109. Access mode: https://link.aps.org/doi/10.1103/PhysRevA.103.053109.

  12. Cox D., Little J., O’Shea D. Ideals, varieties, and algorithms. 3 ed. Springer, 2007.

  13. Zharkov A.Yu., Blinkov Yu.A. Involution approach to solving systems of algebraic equations // Proceed. of the 1993 Inter. IMACS Symp. on Symbolic Comput. Laboratoire d’Informatique Fondamentale de Lille, France, 1993. P. 11–16.

  14. Zharkov A.Yu., Blinkov Yu.A. Solving zero-dimensional involutive systems // Progress in Math. Basel: Birkhauser, 1996. V. 143. P. 389–399.

  15. Blinkov Yu.A. Division and algorithms in the ideal membership problem [Deleniye i algoritmy v zadache o prinadlezhnosti k idealu] // Izvestija Saratovskogo universiteta. 2001. V. 1. № 2. P. 156–167. In Russ.

  16. Apel J. A Gröbner approach to involutive bases // J. of Symbolic Comput. 1995. V. 19. № 5. P. 441–458.

  17. McCarthy J. Recursive functions of symbolic expressions and their computation by machine, Part I // Communicat. of ACM. 1960. V. 3. № 4. P. 184–195.

  18. Jones R., Hosking A., Moss J., Eliot B. The garbage collection handbook: the a of automatic memory management. CRC Applied Algorithms and Data Structures Ser. Chapman and Hall / CRC Press / Taylor & Francis Ltd., 2011. ISBN: 978-1-4200-8279-1.

  19. Клейн Ф. Лекции об икосаэдре и решении уравнений пятой степени. M.: Наука, 1989.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал