ЛОРАНОВЫ РЕШЕНИЯ ЛИНЕЙНЫХ ДИФФЕРЕНЦИАЛЬНЫХ СИСТЕМ
(расширяемое эссе)
С.А.Абрамов, А.А.Рябенко, Д.Е.Хмельнов
$\newcommand{\val}{\mathrm{val} \,}
\newcommand{\zs}{\mathbb{Z}}
\newcommand{\qs}{\mathbb{Q}}
\newcommand{\Mat}{\mathrm{Mat}}
\newcommand{\ddx}{\frac{d}{dx}}
\newcommand{\ind}{\!^{^{\bf ⓘ} }\!}
\newcommand{\ord}{\mathrm{ord} \,}
\newcommand{\rps}{\mbox{“редукция + сдвиг”}}
\newcommand{\Cnstrnts}{\mathcal{C}}
$
Предлагаемое эссе допускает расширения, согласованные с их активными (“кликабельными”) значками на полях текста:
$\,$сказать больше,
$\,$привести пример,
$\,$указать, где об этом написано подробнее,
$\,$задать контрольный вопрос с проверкой ответа,
$\,$закрыть ранее открытое расширение (но такая возможность блокируется, если имеется ссылка на это расширение из другого открытого расширения; тогда значок будет автоматически заменен на ).
Расширяемое эссе первоначально выглядит как отображенный на экране компактный текст, который дает первое представление о предмете. Если содержание какого-то абзаца текста привлекает интерес и если этот абзац помечен на полях значком расширения того или иного типа, то такой абзац можно расширить. Расширенный текст сам может допускать расширения и т.д. От любого выполненного расширения можно отказаться и вернуться к предыдущему варианту текста. Итоговый вариант можно сохранить, распечатать.
Предлагаемое далее эссе тематически связано с вопросами компьютерной алгебры, рассматриваемыми в полугодовом спецкурсе для студентов 2-го года магистратуры факультета ВМК МГУ.
Наряду с “расширяемым эссе”, возможно, уместен термин “расширяемый конспект”.
0.$\,$Введение. Предварительные сведения. Алгоритмы компьютерной алгебры предназначены в основном для компьютерного решения задач с имеющими вид формул исходными данными и результатами.
Ниже обсуждаются начальные аспекты построения решений систем линейных обыкновенных дифференциальных уравнений. Классическая матричная (линейная) алгебра дает удобные средства работы с системами линейных алгебраических уравнений. Мы касаемся вопросов алгебры матриц, элементами которых служат скалярные линейные обыкновенные дифференциальные операторы вида $a_r(x)\frac {d^r}{dx^r}+\dots +a_1(x)\frac {d}{dx}+a_0(x).$ Характер коэффициентов этих скалярных операторов оговаривается специально. Такие матрицы являются удобным средством работы с рассматриваемыми системами.
Далее $K$ — это поле характеристики 0, например (но не обязательно), — некоторое числовое поле. Считаем, что
–$\,$имеются алгоритмы сложения и умножения элементов из $K$ и алгоритм проверки равенства нулю элемента из $K$;
–$\,$кольцо целых чисел вложено в $K$: $\zs\subset K$;
–$\,$имеется алгоритм нахождения целочисленных корней алгебраических уравнений над $K$ с одной неизвестной $n.$
Мы используем для кольца полиномов над $K$ от некоторой переменной $x$ обычное обозначение $K[x],$ для кольца формальных степенных рядов — обозначение $K[ [x] ],$ для поля рациональных функций — обозначение $K(x),$ для поля формальных лорановых рядов — обозначение $K((x)).$
Кольцо полиномов от нескольких переменных $x_1,\dots,x_n$ обозначается через $K[x_1,\dots ,x_n].$
Для кольца квадратных матриц размера $m\times m$ с элементами в некотором кольце или поле $P$ используется обозначение ${\rm Mat}_m(P).$ Для матрицы $M,$ не обязательно квадратной, $M^T$ обозначает транспонированную матрицу. Для квадратной матрицы $M$ под $\det M$ понимается ее определитель.
Множество (кольцо относительно поэлементного сложения и умножения) двусторонних бесконечных последовательностей будет обозначаться через $K^{\infty}.$ Для элементов таких последовательностей будем прибегать к скобочной записи вместо индексной. Например, для элементов последовательности $c\in K^{\infty}$ будем писать $c(n)$ вместо $c_n.$ Нижние индексы используются для элементов конечных наборов — векторов, матриц и т.д.
Рассматривая в этом эссе бесконечные ряды, мы не затрагиваем вопросов сходимости. Все обсуждаемые ряды — формальные. Формальным лорановым рядом над полем $K$ называется выражение (формальная сумма, “картинка”) вида $$\sum _{n=-\infty}^{\infty}c(n)x^n, \tag{1} \label{fsu-}$$ где все значения $c(n)$ (коэффициенты ряда) принадлежат полю $K,$ и при этом для каждого ненулевого формального лоранового ряда $s(x)$ существует принадлежащее множеству целых чисел $\zs$ значение $\nu$ такое, что $c(n)=0$ при $n<\nu.$ Наибольшее такое $\nu$ называется валюацией ряда, для этого числа используется обозначение $\val s(x).$ Для нулевого ряда $s(x)$ по определению $\val s(x)=\infty.$ Ряд $(\ref{fsu-}),$ имеющий валюацию $\nu,$ можно записать как $\sum _{n=\nu}^{\infty}c(n)x^n,$ или $c(\nu)x^{\nu}+c(\nu+1)x^{\nu+1}+\dots$
Формальные лорановы ряды можно складывать и умножать. Эти операции коммутативны и ассоциативны, при этом умножение дистрибутивно относительно сложения. Обычный полином от $x$ может рассматриваться как лоранов ряд (“недостающие” коэффициенты являются нулевыми). Мы можем, например, умножить полином на лоранов ряд.
Там, где это не вызывает недоразумений, мы будем вместо “формальный лоранов ряд” писать “лоранов ряд” или просто “ряд”. Ряд, все коэффициенты которого равны нулю, обозначается символом 0. Выделяют также ряд, все коэффициенты $c(n)$ которого, кроме $c(0),$ равны нулю, при этом $c(0)=1;$ он обозначается символом 1.
Ряды над полем $K$ сами образуют поле $K((x)).$ В поле всех лорановых рядов можно выделить кольцо $K[ [x] ]$ рядов, валюации которых неотрицательны. Это кольцо называется кольцом формальных степенных рядов над $K.$
Производная ряда $s(x)=\sum _{n=-\infty}^{\infty}c(n)x^n$ определяется как $\frac d{dx}s(x)=$ $s'(x)=$ $\sum _{n=-\infty }^{\infty}d(n)x^n,$ где $d(n)=(n+1)c(n+1)$ для всех $n$ (из этого определения следует, что коэффициент $d(-1)$ производной всегда равен нулю). Имеют место обычные соотношения, например, $$(as(x)+bt(x))'=as'(x)+bt'(x),\;\;(s(x)t(x))'=s'(x)t(x)+s(x)t'(x),$$ $s(x),t(x)\in K((x)),\; a,b\in K.$
В дальнейшем нам будет удобным иметь также дело с рядами, коэффициентами которых служат вектор-столбцы: $$y(x)=\sum_{n=-\infty}^{\infty}c(n)x^n, \tag{2} \label{fsu}$$ $c(n)=(c_1(n), \dots, c_m(n))^T.$ При каждом $n\in \zs$ имеем $c(n)\in K^m;$ в этом отличие от $(\ref{fsu-}),$ где $c(n)$ — последовательность с элементами из $K.$ Будем также говорить, что $(\ref{fsu})$ — это вектор-столбец с компонентами в виде рядов: $y(x)=(y_1(x),\dots, y_m(x))^T.$
Решения уравнений и систем, имеющие вид таких лорановых рядов, как $(\ref{fsu-})$ и $(\ref{fsu})$ соответственно, называются лорановыми решениями.
1.$\,$Индуцированные системы. Обсудим ситуацию, когда дифференциальная система обладает лорановым решением в виде двустороннего ряда $(\ref{fsu}).$ Будет показано, что для системы линейных обыкновенных дифференциальных уравнений $$A_r(x)\,y^{(r)}(x)+A_{r-1}(x)\,y^{(r-1)}(x)+\cdots +A_0(x)\,y(x) = 0, \tag{3} \label{i2dfr}$$ где $A_0(x),A_1(x),\dots,A_r(x)\in \Mat_m(K[x]),$ это возможно, если и только если последовательность вектор-столбцов $c(n)=(c_1(n), \dots, c_m(n))^T,$ $n\in \zs,$ удовлетворяет так называемой индуцированной рекуррентной системе $$B_l(n)c(n+l)+B_{l-1}(n)c(n+l-1)+\cdots +B_t(n)c(n+t)=0, \tag{4} \label{indic}$$ где $l,t\in \zs,$ $l\geqslant t$ и $B_l(n),\dots ,B_t(n)\in \Mat_m(K[n]).$ В $(\ref{indic})$ даже не предполагается существование $\nu\in \zs$ такого, что $c(n)=0$ при $n<\nu,$ соотношения $(\ref{indic})$ будут выполняться и при рассмотрении ряда $(\ref{fsu}),$ имеющего ненулевые коэффициенты при $x^n$ с отрицательными $n$ со сколь угодно большой абсолютной величиной.
$A_r(x),$ $B_l(n)$ — ненулевые матрицы, это ведущие матрицы систем $(\ref{i2dfr}),$ $(\ref{indic}).$
Для начала рассмотрим подробности случая одной неизвестной функции, иначе говоря — скалярного случая $m=1.$ Преобразование (замена) $$ x\rightarrow \sigma^{-1},\;\; \ddx\rightarrow (n+1)\sigma \tag{5} \label{dtarr} $$ переводит дифференциальный оператор $$L=a_r(x)\frac {d^r}{dx^r}+\dots +a_1(x)\ddx+a_0(x)$$ в индуцированный рекуррентный оператор $$ L\ind=b_l(n) \sigma^l +\dots +b_t(n)\sigma^t, $$ где $\sigma$ — операция сдвига: $\sigma c(n)$ — это такая последовательность, назовем ее $d(n),$ что $d(n)=c(n+1)$ для всех целых $n;$ соответственно $\sigma ^{-1}$ — сдвиг в обратную сторону: если $d(n)= \sigma^{-1}c(n),$ то $d(n)=c(n-1)$ для всех целых $n.$
При любой двусторонней последовательности $c(n)\in K^{\infty}$ применение $L$ к ряду $(\ref{fsu-})$ даст ряд $\sum _{n=-\infty}^\infty \bar{c}(n)x^n.$ Двусторонняя последовательность $\bar{c}(n)$ получается применением $L\ind$ к последовательности $c(n)\!:$ $$\bar{c}(n)=b_l(n)c(n+l)+b_{l-1}(n)c(n+l-1)+\dots+b_t(n)c(n+t), \;\; n\in \zs.$$
Мы имеем, таким образом, $$L(y)=0\; \Longleftrightarrow \;L\ind \,(c)=0.$$ Это применимо и к общему случаю $m\geqslant 1.$
Обратимся к иллюстративному примеру.
Пусть, например, $L=x\frac d{dx}-(x-1),$ тогда $L\ind=\sigma^{-1}(n+1)\sigma-(\sigma^{-1}-1)=$ $n+1-\sigma^{-1}.$ Если мы интересуемся лорановыми решениями уравнения $L(y)=0$ с указанным $L,$ то вид $L\ind$ позволяет заметить, что $c(n)=0$ при $n\leqslant -2,$ а $c(-1)$ может быть выбрано произвольно, и $c(n)=\frac {c(-1)}{(n+1)!}$ при $n\geqslant -1.$
Для множества скалярных дифференциальных операторов мы используем обозначение $K[x, \frac d{dx}],$ так как эти операторы выглядят как полиномы над $K$ от $x$ и $\frac d{dx}.$ Эти “полиномы” считаются записанными по степеням $\frac d{dx}.$ Аналогично для случая рекуррентных операторов используется обозначение $K[n, \sigma, \sigma^{-1}].$
Общая теория такого рода “некоммутативных полиномов” была предложена О.$\,$Оре.
В общем случае исходную систему можно записать в матричной форме $$\left(\begin{array}{ccc} L_{11}&\dots & L_{1m} \cr \dots&\dots & \dots \cr L_{m1}&\dots & L_{mm} \end{array}\right)y(x)=0, \tag{6} \label{om}$$ $L_{ij}\in K[x,\frac d{dx}],$ $i,j=1,\dots,m.$
Индуцированная рекуррентная система для $(\ref{om})$ имеет вид $$\left(\begin{matrix} L\ind_{11}&\dots & L\ind_{1m}\cr \dots&\dots & \dots\cr L\ind_{m1}&\dots & L\ind_{mm} \end{matrix}\right) c(n)=0, \tag{7} \label{om+}$$ $L\ind_{ij}\in K[n,\sigma,\sigma^{-1}],$ $i,j=1,\dots,m.$ Каждый оператор $L\ind_{ij}$ получается из $L_{ij}$ преобразованием $(\ref{dtarr}).$
2.$\,$Операторные матрицы. Система $(\ref{i2dfr})$ может быть записана как $L(y)=0,$ где $L$ — дифференциальный оператор с матричными коэффициентами: $$L = A_r(x)\frac {d^r} {dx^r}+A_{r-1}(x)\frac {d^{r-1}} {dx^{r-1}}+\cdots+A_0(x).\tag{8} \label{i2op}$$ Понятие ведущей матрицы оператора определяется так же, как для соответствующей системы; $r$ является порядком оператора $L$ (пишем $r=\ord L$).
При $m=1$ оператор $L$ становится скалярным оператором.
Оператор $(\ref{i2op})$ может быть представлен единой матрицей, принадлежащей $\Mat _m(K[x, \ddx])$ (см. $(\ref{om})$): $$L = \left( \begin{array}{lll} L_{11}&\dots&L_{1m}\cr \dots&\dots&\dots\cr L_{m1}&\dots&L_{mm}\cr \end{array} \right)\!\!, \tag{9} \label{mrep}$$ $ L_{ij}\in K [x, \ddx],$ $i,j=1,\dots,m,$ при этом $\max_{i,j}\ord L_{ij}=r.$
Таким образом, мы имеем представления $(\ref{i2op})$ и $(\ref{mrep})$ рассматриваемых операторов. Будем выбирать форму представления, исходя из соображений удобства.
Матрицы, не обязательно квадратные, с принадлежащими $K[x,\ddx]$ элементами, мы будем называть операторными матрицами. В тех случаях, когда это не вызывает недоразумений, будем использовать термин оператор для принадлежащих $\Mat _m(K[x, \ddx])$ квадратных операторных матриц, а термин “матрица” — для матриц с элементами из $K[x].$
Матрица $A_r(x)$ является ведущей матрицей системы $L(y)=0$ и оператора $L$ при любой выбираемой форме представления системы и оператора. Аналогичные определения принимаются для рекуррентных систем и операторов. Для рекуррентной системы $L\ind(c)=0$ с оператором $L\ind$ вида $$B_l(n)\sigma^{l} + B_{l-1}(n)\sigma^{l-1}+ \cdots +B_t(n)\sigma^{t} \tag{10} \label{i2sp} $$ вводятся понятия ведущей $B_l(n) \neq0$ и трейлинговой $B_t(n)\neq 0$ матриц, ведущего $l$ и трейлингового $t$ порядков. Порядок рекуррентного оператора: $\ord L\ind = l-t.$
3.$\,$Независимость уравнений. Для $i$-й строки произвольной матрицы $M$ будем использовать обозначение $M_{i,*}.$ Если исходная система $L(y)=0$ записана в виде $(\ref{om}),$ то условие независимости над $K[x,\frac d{dx}]$ уравнений исходной системы равносильно независимости строк операторной матрицы $L.$ Определив результат умножения скалярного оператора из $K[x,\frac d{dx}]$ на строку операторной матрицы как строку, составленную из произведений этого оператора на соответствующие элементы строки, мы можем записать условие независимости уравнений системы в форме условия независимости строк $L$ над $K[x,\frac d{dx}]\!:$ если $F_1L_{1,*}+\dots+F_mL_{m,*}=0,$ для некоторых $F_1,\dots,F_m\in K[x,\frac d{dx}],$ то $F_1=\dots=F_m=0.$
Операторная матрица (оператор) имеет полный ранг, если все строки этой матрицы независимы в указанном смысле.
Подобным образом вводится понятие независимости строк над некоторым кольцом $R$ для любой матрицы $M$ с элементами из $R.$
Преобразование $(\ref{dtarr})$ задает гомоморфизм кольца дифференциальных операторов $K[x,\frac d{dx}]$ в кольцо $K[n,\sigma,\sigma^{-1}]$ рекуррентных операторов. Будем использовать для этого гомоморфизма обозначение ${\cal M}_{\frac d{dx}}.$ Формулы $(\ref{dtarr})$ можно переписать как $${\cal M}_{\frac d{dx}}(x)=\sigma^{-1},\;\; {\cal M}_{\frac d{dx}}\left(\frac d{dx}\right)=(n+1)\sigma. $$ Этот гомоморфизм расширяется на кольцо $K[x, x^{-1},\frac d{dx}]\!:$ $${\cal M}_{\frac d{dx}}(x^{-1})=\sigma,$$ при этом он становится изоморфизмом $${\cal M}_{\frac d{dx}} \,: \, K[x, x^{-1},\frac d{dx}] \rightarrow K [n,\sigma,\sigma^{-1}]. \tag{11} \label{isom}$$ Обратный изоморфизм определяется формулами
$${\cal M}_{\frac d{dx}}^{-1}(n)= x\frac d{dx},\;\;{\cal M}_{\frac d{dx}}^{-1}(\sigma)= x^{-1},\;\; {\cal M}_{\frac d{dx}}^{-1}(\sigma^{-1})=x.$$
Теорема$\,$1. Пусть уравнения дифференциальной системы $L(y)=0$ независимы над $K[x,\frac d{dx}],$ тогда уравнения индуцированной рекуррентной системы $L\ind (c)=0,$ $L\ind={\cal M}_{\frac d{dx}}L,$ независимы над $K[n,\sigma, \sigma^{-1}].$
Далее мы увидим, что существенным удобством при нахождении решений дифференциальных систем служит невырожденность ведущей матрицы оператора $L\ind,$ т.е. матрицы $B_l(n)$ в $(\ref{indic}).$
4.$\,$Определяющие уравнения. Напомним, что для ненулевого степенного или лоранова ряда $s(x)=\sum_n c(n)x^n$ его валюацией, обозначаемой через $\val s(x),$ называется целое число $$\min \{n\in \zs \mid c(n)\neq 0\},$$ при этом $\val s(x)=\infty$ для нулевого ряда $s(x).$ Валюацией вектора или матрицы с элементами-рядами называется наименьшая из валюаций элементов этого вектора или матрицы.
В классической теории обыкновенных дифференциальных уравнений для нахождения валюаций (порядков нулей и полюсов) возможных аналитических решений в скалярном случае применяются алгебраические уравнения, называемые определяющими.
Для нахождения нижних границ валюаций принадлежащих $K((x))^m$ решений дифференциальных систем нам нужны какие-то варианты определяющих уравнений. Индуцированные рекуррентные системы позволяют построить такого рода уравнение.
Пусть $(\ref{indic})$ является индуцированной системой для $(\ref{i2dfr}).$ Если для некоторого целого $n_0$ выполнено $\det B_l(n_0)\neq 0,$ то домножением $(\ref{indic})$ на $(B_l(n_0))^{-1}$ получаем $$ \begin{split}c(n_0+l)=-(B_l(n_0))^{-1}(B_{l-1}(n_0)c(n_0+l-1) & +\dots  + \\[0.3em] &\; +\, B_t(n_0)c(n_0+t)), \end{split} \tag{12} \label{vyr}$$ и при известных значениях $c(n_0+l-1),\dots, c(n_0+t)$ значение $c(n_0+l)$ однозначно определяется. В частности, при $c(n_0+l-1)=\dots =c(n_0+t)=0$ имеем $c(n_0+l)=0.$ Отсюда, если $\det B_l(n)$ — ненулевой полином от $n,$ то для валюации $\nu$ какого-либо лоранова решения системы $(\ref{i2dfr})$ должно выполняться $\det B_l(\nu -l)=0.$
Итак, при $\det B_l(n_0)\neq 0$ значение $c(n_0+l)$ однозначно определяется по предшествующим элементам последовательности $c(n)$ вектор-столбцов. Если же $\det B_l(n_0)= 0,$ то вектор-столбец $c(n_0+l)$ удовлетворяет системе линейных алгебраических уравнений $$B_l(n_0)c(n_0+l)=-(B_{l-1}(n_0)c(n_0+l-1)+\dots +B_t(n_0)c(n_0+t))$$ с вырожденной матрицей $B_l(n_0).$ Когда эта система неоднородна, то она в каких-то случаях имеет решение в $K^m,$ а в каких-то — нет.
Получаем, в частности, следующее
Предложение$\,$1. Пусть $y(x)=(y_1(x),\dots ,y_m(x))^T \in K((x))^m$ — решение дифференциальной системы, для которой $(\ref{indic})$ является индуцированной рекуррентной системой. Тогда для значения $\nu=\val y(x)$ выполняется $$\det {B}_l(\nu -l)=0.\tag{13}\label{ind2_1}$$
Сразу заметим, что если матрица $B_l(n)$ вырождена, т.е. ее определитель — это нулевой полином от $n,$ то это предложение ничего не дает для определения возможных значений $\val y(x),$ так как уравнение $(\ref{ind2_1})$ превращается в бесполезное тождество $0=0.$ Выходом могло бы быть приведение системы $L\ind(c)=0$ к виду, в котором ведущая матрица невырождена и, вместе с этим, либо преобразованная система имеет решения, совпадающие с решениями системы $L\ind(c)=0,$ либо, по крайней мере, ни одно решение исходной системы не будет потеряно, хотя появление лишних решений возможно.
Ниже мы дадим алгоритм преобразования заданной рекуррентной системы $R(c)=0$ с оператором $R\in K[n,\sigma,\sigma^{-1}]$ в систему ${\bar R}(c)=0$ с невырожденной ведущей матрицей, это преобразование не вызывает потери каких-либо решений первоначальной системы. Получаемая в результате система будет называться охватывающей системой для системы $R(c)=0,$ а оператор ${\bar R}$ — охватывающим оператором для $R.$
С помощью охватывающей для $(\ref{indic})$ системы $$ \bar{B}_l(n)c(n+l)+ \bar{B}_{l-1}(n)c(n+l-1)+\cdots +\bar{B}_t(n)c(n+t)=0\tag{14} \label{lindic}$$ и следующей теоремы можно будет найти некоторую нижнюю границу валюаций лорановых решений исходной дифференциальной системы $(\ref{i2dfr}).$
Теорема$\,$2. Пусть $y(x)=(y_1(x),\dots ,y_m(x))^T \in K((x))^m$ — решение дифференциальной системы, для которой $(\ref{indic})$ является индуцированной рекуррентной системой, а $(\ref{lindic})$ — охватывающей для $(\ref{indic}).$ Тогда для значения $\nu=\val y(x)$ выполняется $$\det \bar{B}_l(\nu -l)=0.\tag{15}\label{ind2}$$
Условимся называть уравнения вида $(\ref{ind2})$ для $\nu$ определяющими уравнениями исходной дифференциальной системы.
Решение рекуррентной системы, имеющее вид последовательности, далее будем называть секвенциальным решением.
5.$\,$Алгоритм EG$_{\sigma}.$ Алгоритм EG$_{\sigma}$ позволяет получить для данного оператора $R\in \Mat _m(K[n, \sigma, \sigma^{-1} ])$ полного ранга охватывающий оператор $\bar{R}=MR,$ имеющий невырожденную ведущую матрицу, при этом $M$ — некоторый унимодулярный (т.е. имеющий обратный в $\Mat _m(K(n)[\sigma, \sigma^{-1} ])$) оператор; построение оператора $M$ не выполняется, но при необходимости оно может быть включено в алгоритм.
Пусть оператор полного ранга $R$ имеет вид $(\ref{i2sp}).$ Если $1\leqslant i\leqslant m,$ то определим ведущий порядок $\alpha _i(R)$ строки с номером $i$ как наибольшее целое $k,$ $l\geqslant k\geqslant t,$ при котором $i$-я строка матрицы $B_k(n),$ т.е. строка $(B_k)_{i,*},$ ненулевая; соответственно, трейлинговым порядком $\beta _i(R)$ этой строки назовем наименьшее целое $k,$ $l\geqslant k\geqslant t,$ при котором $i$-я строка матрицы $B_k(n),$ т.е. строка $(B_k)_{i,*},$ ненулевая. Когда это не вызывает недоразумений, будем вместо $\alpha _i(R),$ $\beta _i(R)$ писать $\alpha _i,$ $\beta _i.$
Оператор $(\ref{i2sp})$ можно представить как совокупность операторов вида $$R_i=(B_l(n))_{i,*}\ \sigma ^l+(B_{l-1}(n))_{i,*}\ \sigma ^{l-1}+\dots +(B_t(n))_{i,*}\ \sigma ^t, \tag{16} \label{rr} $$ $i=1,\dots, m.$ Видно, что $R_i(c)$ дает $i$-ю компоненту вектор-столбца $R(c).$
И теперь — сам алгоритм:
Проверить, являются ли строки ведущей матрицы оператора $R$ линейно зависимыми над $K[n].$ Если нет, то оператор $R$ не изменяется и алгоритм останавливается. В противном случае следует серия шагов, каждый из которых состоит из двух этапов.
Этап редукции: Вычислить коэффициенты $p_1(n), \dots, p_m(n) \in K[n]$ зависимости строк ведущей матрицы. Взять какое-нибудь $i,$ $1\leqslant i\leqslant m,$ для которого $$\beta_i=\min_{ {1\leqslant j\leqslant m, p_j\neq 0} }\; \beta_j.$$ Выполнить $R_i:=\sum _{k=1}^m p_k(n)R_k.$ В ведущей матрице оператора $R$ при этом $i$-я строка становится нулевой, т.е. $(B_l(n))_{i,*}$ в $R_i$ становится нулевой строкой. Имеем неравенство $\alpha _i<l.$
Этап сдвига: К $R_i$ (значение $i$ выбрано на этапе редукции) применить $\sigma^{l-\alpha _i},$ имеем $\alpha _i=l.$
После некоторого конечного числа шагов $\rps$ получаем оператор с невырожденной ведущей матрицей.
Завершимость гарантируется тем, что на каждом шаге $\rps$ этап редукции не уменьшает $\beta _i,$ а этап сдвига увеличивает $\beta _i.$ Границы изменения значения $\sum _{j=1}^m\beta _j$ устанавливаются несложно.
Алгоритм EG$_{\sigma}$ может быть применен не только к оператору $R$ вида $(\ref{i2sp}),$ но и непосредственно к системе $R(c)=0.$ Вместо указанного в $(\ref{rr})$ оператора $R_i,$ можно рассматривать $i$-е уравнение нашей системы, которое запишется как $R_i(c)=0.$
В следующем разделе мы сделаем добавление к алгоритму EG$_{\sigma}.$
6.$\,$Линейные ограничения. Если $R$ имеет вид $(\ref{i2sp}),$ то предлагаемый дополненный вариант алгоритма EG$_{\sigma}$ позволяет получить для системы $R(c)=0$ вместе с охватывающей системой, еще и конечное множество соотношений, называемых линейными ограничениями, каждое из которых затрагивает конечное число элементов любого из секвенциальных решений исходной системы. Когда для секвенциального решения охватывающей системы удовлетворены все линейные ограничения, это решение является решением исходной рекуррентной системы.
Если на некотором этапе редукции в ходе выполнения алгоритма EG$_{\sigma}$ $i$-е уравнение $$R_{i} (c) = \sum_{j=1}^m \sum_{k=t}^l (B_k(n))_{i,j}\, c_j(n+k) = 0 \tag{17} \label{s174} $$ заменяется в системе линейной комбинацией уравнений с коэффициентами, являющимися полиномами $p_1(n),\dots,p_m(n),$ то изменяемое $i$-е уравнение домножается на полином $p_i(n),$ который может обращаться в $0$ при некоторых значениях $n.$ Каждое такое значение $n=\eta$ подставляется вместо $n$ в уравнение $(\ref{s174}).$ Получаемое линейное ограничение имеет вид линейного соотношения с коэффициентами из $K$ для значений $c_j(\eta+k)$ ($j=1,\dots, m$; $k=t,t+1,\dots,l$): $$\sum_{j=1}^m \sum_{k=t}^l (B_k(\eta))_{i,j}\, c_j(\eta+k) = 0.$$
7.$\,$Полиномиальные коэффициенты: лорановы решения. Для имеющей полиномиальные коэффициенты системы $(\ref{i2dfr})$ обсудим поиск таких ее решений, компоненты которых являются лорановыми рядами: $$ y(x) = \sum_{n=\nu}^\infty c(n)\,x^n, \;\; \nu\in \zs, $$ где $c(n)=(c_1(n),\dots ,c_m(n))^T\in K^m,$ $n \in \zs.$ Как договаривались, эти решения мы называем лорановыми. Вопрос о сходимости рядов мы не рассматриваем.
Если для $(\ref{i2dfr})$ индуцированная система $(\ref{indic})$ имеет вырожденную ведущую матрицу $B_l(n),$ то мы с помощью EG$_{\sigma}$ строим для индуцированной системы охватывающую рекуррентную систему $(\ref{lindic}).$ При этом возникает конечное множество — не исключается, что пустое — линейных ограничений.
Пользуясь теоремой 2, находим нижнюю границу $n_0$ валюаций решений, где $n_0$ — наименьший целый корень определяющего уравнения $\det {\bar {B}}_l(n-l)=0$; если целых корней нет, то система не имеет ненулевых лорановых решений.
Ряд представляется начальными членами $$c(\nu)\, x^{\nu}+\dots+ c(N)\, x^{N}, \tag{18} \label{inte}$$ и значение $N,$ большее или равное $\nu,$ выбирается так, что, во-первых, вычисление последующих членов $c(n)\,x^n$ уже не требует учета линейных ограничений, и, во-вторых, определитель ведущей матрицы охватывающей рекуррентной системы в ходе вычисления этих последующих членов не обращается в нуль.
Коэффициенты начальных членов $(\ref{inte})$ можно найти из системы линейных алгебраических уравнений, построенной с помощью охватывающей индуцированной рекуррентной системы с добавлением линейных ограничений, в которые подставлены нули вместо неизвестных $c_1(n),\dots, c_m(n)$ при $n<\nu.$
Любое количество последующих членов может быть получено с помощью охватывающей рекуррентной системы, т.е. по формуле $$c(n)=-\bar{B}^{-1}_l(n-l)\left(\bar{B}_{l-1}(n-l)c(n-1)+\cdots +\bar{B}_t(n-l)c(n-l+t)\right),$$ $n=N+1, N+2,\dots$
8.$\,$Коэффициенты-ряды: полнота ранга. В качестве поля коэффициентов дифференциальных систем в наших рассмотрениях теперь будет выступать $K((x)).$ По прежнему $K$ имеет характеристику $0.$ В $K((x))$ определено обычное дифференцирование (см. 0. Введение. Предварительные сведения.)
Дифференциальные системы будет удобно записывать с помощью операции $\theta=x\frac d{dx}$ вместо обычной операции дифференцирования $\frac d{dx}$ (переход от одной формы записи к другой выполняется легко): \begin{equation} A_r(x)\theta^{r}y(x)+A_{r-1}(x)\theta^{r-1}y(x)+\cdots +A_0(x)y(x) = 0. \tag{19} \label{i2m} \end{equation}Как и прежде, $y(x)=(y_1(x),y_2(x),\dots, y_m(x))^T$ — это вектор-столбец неизвестных лорановых рядов от $x.$ Относительно матриц \begin{equation}A_0(x),A_1(x),\dots ,A_r(x) \tag{20}\label{10matcoef} \end{equation} считаем, что $A_i(x)\in \Mat_m (K[ [x] ]),$ $i=0,1,\dots,r,$ и $A_r(x)$ — ненулевая ведущая матрица. Число $r$ является порядком системы.
Через ${\cal D} _m$ обозначим кольцо дифференциальных операторов с матричными коэффициентами $\Mat _m(K[ [x] ])[\theta].$ Систему $(\ref{i2m})$ удобно будет записывать как $L(y)=0,$ где \begin{equation} L=A_r(x)\theta^{r}+A_{r-1}(x)\theta^{r-1}+\cdots +A_0(x)\in {\cal D} _m. \tag{21} \label{Loper} \end{equation} Оператор $L$ представляется матрицей $$\left( \begin{array}{lll} L_{11}&\dots&L_{1m}\\ \dots&\dots&\dots\\ L_{m1}&\dots&L_{mm}\\ \end{array} \right)\!, \tag{22} \label{10mrep}$$ в которой $L_{ij}\in K[ [x] ] [\theta],$ $i,j=1,\dots,m,$ при этом $\max_{i,j}\ord L_{ij}=r.$ Система $L(y)=0$ и оператор $L$ имеют полный ранг, если и только если строки матрицы $(\ref{10mrep})$ линейно независимы над $K[ [x] ][\theta]$ (см. 3. Независимость уравнений.).
9.$\,$Представление бесконечных рядов. Бесконечные последовательности коэффициентов рассматриваемых рядов задаются алгоритмически. Пусть система $L(y)=0$ имеет вид $(\ref{i2m})$ и при этом для каждого элемента $a(x)$ какой-либо из матриц указан такой алгоритм $\Xi _a,$ что $a(x)=\sum_{k=0}^{\infty}\Xi_a(k)x^k.$ Согласно классическому результату Тьюринга, мы в общем случае не можем проверить, является ли заданный ряд нулевым или нет.
Мы будем предполагать, что заданная система имеет полный ранг. Требование полноты ранга системы является существенным: отказ от него делает неразрешимой задачу алгоритмической проверки существования решений.
Предложение$\,$2. Задача проверки для произвольной заданной системы вида $(\ref{i2m})$ наличия у нее решения в $K((x))^m\setminus \{0\}$ алгоритмически неразрешима. Алгоритмически неразрешима и задача проверки полноты ранга данной системы.
Как будет показано, если относительно системы вида $(\ref{i2m})$ заранее известно, что она имеет полный ранг, то все ее лорановы решения могут быть найдены алгоритмически.
10.$\,$Бесконечные индуцированные рекуррентные системы. Итак, предполагается, что система имеет полный ранг. Дополнительно будем предполагать, что по крайней мере один из коэффициентов-рядов в каждом уравнении системы имеет ненулевой свободный член, т.е. имеет валюацию $0.$
Система $L(y)=0$ имеет лораново решение \begin{equation}y(x)=c(\nu)x^{\nu}+c(\nu+1)x^{\nu+1}+\dots, \tag{23} \label{LSol}\end{equation} если и только если последовательность ${c(n)=(c_1(n),\dots ,c_m(n))^T\in K^m},$ ${n \in \zs},$ векторных коэффициентов решения $y(x),$ т.е. векторная последовательность \begin{equation}\dots , 0,0,c(\nu),c(\nu+1),\dots, \tag{24} \label{LSol_coeffs}\end{equation} удовлетворяет индуцированной рекуррентной системе. В случае, когда элементы коэффициентов-матриц $(\ref{10matcoef})$ суть ряды, индуцированная система имеет бесконечный порядок (в отличие от уже рассмотренного случая системы с полиномиальными коэффициентами, когда порядок индуцированной системы конечен).
Преобразование \begin{equation}x\rightarrow \sigma ^{-1},\;\;x ^{-1}\rightarrow \sigma ,\;\;\theta \rightarrow n \tag{25} \label{ismr}\end{equation} переводит исходную дифференциальную систему $L(y)=0$ в индуцированную рекуррентную систему \begin{equation}B_0(n)c(n)+ B_{-1}(n)c(n-1)+ \cdots =0, \tag{26} \label{i2i}\end{equation} где
–$\,B_0(n),B_{-1}(n),\dots \in \Mat _m (K[n])$ и каждый из полиномиальных элементов этих матриц имеет степень меньшую или равную $r$;
–$\,B_0(n)$ является ненулевой матрицей, это ведущая матрица системы $(\ref{i2i});$
–$\,c(n)=(c_1(n),\dots ,c_m(n))^T$ — такой вектор-столбец неизвестных последовательностей, что $c_i(n)=0$ для всех отрицательных целых $n$ с достаточно большим значением $|n|,$ $i=1,\dots,m.$
Примечание. Ведущий порядок индуцированной рекуррентной системы $(\ref{i2i})$ равен $0.$ Это является следствием того, что дифференциальная система записана с помощью операции $\theta$ и валюация всей системы равна $0.$ Элементы матрицы $B_{-i}(n)$ индуцированной системы целиком определяются коэффициентами при $x^i$ рядов, являющихся элементами матриц $A_r(x),$ $A_{r-1}(x), \dots,$ $A_0(x)$ в исходной дифференциальной системе $(\ref{i2m}).$ Это играет существенную роль при работе с бесконечными рекуррентными системами.
Преобразование $(\ref{ismr})$ задает изоморфизм колец \begin{equation}{\cal M}_\theta:\ \Mat _m(K((x)))[\theta]\rightarrow \Mat _m(K[n])((\sigma^{-1})) \tag{27} \label{mell}\end{equation} (аналогично изоморфизму $(\ref{isom})$ для операторов с полиномиальными коэффициентами). Это преобразование переводит дифференциальный оператор $L$ в индуцированный рекуррентный оператор $$R = B_0(n)+ B_{-1}(n)\sigma^{-1}+ B_{-2}(n)\sigma^{-2}+\dots \tag{28} \label{ind_op}$$ Кольцо $\Mat _m(K[n])[ [\sigma^{-1}] ]$ рекуррентных операторов будет обозначаться через ${\cal E} _m.$
Очевидно, что исходная дифференциальная система $L(y)=0$ имеет полный ранг, если и только если индуцированная система $(\ref{i2i})$ имеет полный ранг, т.е. уравнения системы $(\ref{i2i})$ независимы над $K[n][ [\sigma ^{-1}] ].$
Итак, система $L(y)=0$ имеет лораново решение $(\ref{LSol}),$ если и только если двусторонняя последовательность $(\ref{LSol_coeffs})$ удовлетворяет индуцированной системе $(\ref{i2i}),$ т.е. выполняется \begin{eqnarray*} & &B_0(\nu)c(\nu)=0,\\ & &B_0(\nu+1)c(\nu+1)+B_{-1}(\nu+1)c(\nu)=0,\\ & &B_0(\nu+2)c(\nu+2)+B_{-1}(\nu+2)c(\nu+1)+B_{-2}(\nu+2)c(\nu)=0,\\ & &\dots \end{eqnarray*} Если $\det B_0(n)\neq 0,$ то этот определитель может рассматриваться как определяющее уравнение исходной системы $L(y)=0$: множество целых корней алгебраического уравнения $\det B_0(n)=0$ конечно и содержит множество всех возможных валюаций $\nu$ лорановых решений системы $L(y)=0.$ Однако во многих случаях матрица $B_0(n)$ оказывается вырожденной, даже когда ведущая матрица $A_r(x)$ системы $L(y)=0$ обратима в $\Mat_m(K((x))).$
В разд. 5 обсуждался алгоритм EG$_{\sigma},$ с его помощью после некоторого конечного числа шагов $\rps$ исходная индуцированная система $R(c)=0$ конечного порядка преобразуется в охватывающую систему ${\bar R}(c)=0$: ее ведущая матрица невырождена и множество решений исходной индуцированной системы является подмножеством всех решений охватывающей системы.
Согласно алгоритму EG$_{\sigma},$ исключения (редукции) проводятся в строках с наименьшими трейлинговыми порядками. Но если, как в $(\ref{i2i}),$ трейлинговый порядок равен $-\infty,$ то у нас, вообще говоря, нет возможности выбирать для редукции строку с наименьшим трейлинговым порядком. Требуется некоторая модификация как самого алгоритма, так и доказательства его завершимости.
11.$\,$Алгоритм EG$_{\sigma}^{\infty}.$ Предложим некоторую модификацию последовательности шагов $\rps$ алгоритма EG$_{\sigma}.$ Вместе с преобразованием индуцированной системы мы будем изменять вектор $\gamma=(\gamma_1,\gamma_2,\dots,\gamma_m)$ с положительными целыми элементами. Первоначально $\gamma=(r,\dots,r).$
Шаг $\rps$ выглядит так:
Этап редукции: Использовать любой имеющийся алгоритм, чтобы проверить, зависимы ли строки ведущей матрицы линейно над $K(n),$ и если зависимы, то найти коэффициенты $$p_1(n),\dots ,p_m(n)\in K[n] \tag{29} \label{depnd}$$ зависимости. Положить $$\mu=\max_{{1\leqslant j\leqslant m},\; {p_j(n)\neq 0}}\;\; \{\gamma _j+\deg p_j(n)\}. \tag{30} \label{defmu} $$ Взять $i$ таким, что $$1\leqslant i\leqslant m,\;\;p_i(n)\neq 0,\;\; \gamma _i+\deg p_i(n)=\mu. \tag{31} \label{gamom}$$
Заменить $i$-е уравнение индуцированной рекуррентной системы линейной комбинацией всех ее уравнений с коэффициентами $p_1(n),\dots ,p_m(n).$ В результате, $i$-я строка ведущей матрицы становится нулевой.
Этап сдвига: Применить $\sigma$ к $i$-му уравнению системы.
Увеличить $\gamma_i$ на $\deg p_i(n),$ т.е. положить $\gamma_i:=\mu$ (см. $(\ref{defmu})).$
Последовательность шагов $\rps$ продолжается, пока строки ведущей матрицы остаются линейно зависимыми. Присутствие в алгоритме вектора $\gamma$ существенно используется в доказательстве завершимости. При этом мы никогда не получим нулевого уравнения в рекуррентной системе, потому что уравнения исходной системы линейно независимы над $K[n][ [\sigma ^{-1}_n] ].$ Они остаются линейно независимыми после каждого из шагов $\rps.$
Очевидно, что этап сдвига $i$-го уравнения является эквивалентным преобразованием рекуррентной системы.
Но если полином $p_i(n)$ в $(\ref{depnd})$ имеет целые корни, то этап редукции $i$-го уравнения может породить лишние решения. В этом случае эквивалентность достигается добавлением к полученной рекуррентной системе конечного множества линейных ограничений.
На этапе редукции можно получать множество линейных ограничений, которое мы обозначим посредством $\Cnstrnts.$ Предположим, что мы заменяем $i$-е уравнение системы линейной комбинацией всех уравнений с коэффициентами $(\ref{depnd}),$ при этом $\eta$ является целым корнем полинома $p_i(n).$ Пусть перед этапом редукции $i$-е уравнение имело вид: $$R_i(c) = (B_0(n))_{i,*} \ c(n)+(B_{-1}(n))_{i,*} \ c(n-1)+\dots=0.$$ Если $y(x)=\sum _{n=\nu}^{\infty}c(n)x^n$ служит решением исходной дифференциальной системы, то мы получаем линейное ограничение \begin{equation}\label{cnstr}(B_0(\eta))_{i,*} \ c(\eta)+(B_{-1}(\eta))_{i,*} \ c(\eta-1)+\dots +(B_{-\eta+\nu}(\eta))_{i,*} \ c(\nu)=0.\end{equation}
Здесь существенно, что множество линейных ограничений конечно, и каждое из этих ограничений содержит лишь конечное число ненулевых слагаемых (благодаря тому, что нас интересуют только решения $c(n),$ для которых $c_i(n)=0$ при всех $n,$ меньших нижней границы валюаций лорановых решений исходной дифференциальной системы).
12.$\,$Ленивые вычисления. В основании алгоритма EG$_{\sigma}^{\infty}$ лежит процедура выполнения шагов $\rps.$ Но система $(\ref{i2i}),$ которая должна быть преобразована этой процедурой, имеет бесконечный порядок. Выручают ленивые вычисления с хранением RS-последовательности, содержащей необходимые параметры ранее выполненных редукций и сдвигов.
Для целого $t\geqslant 0$ назовем $t$-усечением системы $(\ref{i2i})$ систему с трейлинговым порядком $t,$ все $t+1$ матричные коэффициенты которой совпадают с начальными $t+1$ матричными коэффициентами $(\ref{i2i})$: $$B_0(n)c(n) +B_{-1}(n)c(n-1) +\dots + B_{-t}(n) c(n-t)=0.$$
Вычисления начинаются с пустой RS-последовательности и $t$-усечения исходной системы, где $t=1.$
Если ведущая матрица $t$-усечения вырождена, для усечения выполняется шаг $\rps.$ Номер $i$ и коэффициенты $(\ref{depnd})$ сохраняются в RS-последовательности. Если оказывается, что на этапе сдвига для вычисления $i$-й строки ведущей матрицы необходимы элементы от $(t+1)$-го матричного коэффициента исходной индуцированной системы (т.е. нам не хватило взятого усечения индуцированной системы), то мы удлиняем усечение, положив $t:= t+1,$ и в получившейся рекуррентной системе повторяем ранее выполненные шаги $\rps$ из RS-последовательности. После этого выполняем шаг сдвига, который потребовал удлинения усечения индуцированной системы и продолжаем дальнейшее выполнение шагов $\rps,$ дополняя RS-последовательность и при необходимости аналогично выполняя дополнительные удлинения этого усечения. Поскольку процесс $\rps$ конечен, то описанные вычисления завершатся. Будет достаточно вычислить конечное число матричных коэффициентов исходной бесконечной индуцированной системы, чтобы получить преобразованную невырожденную ведущую матрицу.
Используя такой же подход, мы можем вычислить необходимое количество следующих матричных коэффициентов преобразованной индуцированной системы, определив по конечной RS-последовательности, какое усечение исходной индуцированной системы для этого будет достаточно, и повторив шаги $\rps$ из RS-последовательности. Также вычисляются и линейные ограничения. Целые корни полиномов из линейных зависимостей используемых в редукциях (эти полиномы сохраняются в RS-последовательности) определяют минимально необходимый размер усечения для вычисления линейных ограничений.
13.$\,$Реализация в системе компьютерной алгебры Maple. Алгоритмы построения решений для систем линейных обыкновенных дифференциальных уравнений реализованы в системе компьютерной алгебры Maple.
На рисунке ниже приведена сессия в системе Maple 2023 построения для системы с полиномиальными коэффициентами индуцированной системы, охватывающей для индуцированной и базиса пространства лорановых решений.
Пакет LFS находится в свободном доступе на веб-станице http://www.ccas.ru/ca/lfs. Там же можно найти Maple-сессии с примерами использования процедур пакета и ссылки на статьи с описанием интерфейса.