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

1. ВВЕДЕНИЕ

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

Среди спектральных методов выделяются два, которые обладают уникальными свойствами, выделяющими их из всех прочих. Это разложение периодических функций в ряды Фурье и разложение функций, аналитических в некоторой окрестности отрезка, по (смещенным) полиномам Чебышёва (что можно рассматривать как модификацию метода Фурье).

Принципиальное отличие этих двух методов от всех прочих спектральных состоит в том, что они одновременно являются также и сеточными, что гарантируется дискретным преобразованием Фурье или Фурье–Чебышёва.

Эти преобразования, которые также можно записать в матричной форме в явном виде, всегда хорошо обусловлены, являются (почти) инволютивными, и не добавляют погрешности к результату при вычислениях с конечной разрядной сеткой (см. [1]).

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

Эти свойства позволяют по таблице полученного решения сделать заключение об аналитических свойствах решения, которые редко заранее известны. Иными словами, скорость убывания коэффициентов Фурье–(Чебышёва) находится в хорошо известной зависимости от гладкости решения. Например, для аналитических функций коэффициенты Фурье–(Чебышёва) убывают экспоненциально быстро.

Для решений, сингулярных в концах интервала, коэффициенты Фурье–(Чебышёва) убывают обычно все же достаточно быстро, чтобы метод был эффективным.

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

Изложенные соображения были положены К.И. Бабенко в основу т.н. методов без насыщения (см. [2]).

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

Применение спектральных методов часто сопряжено с большими техническими трудностями, связанными с учетом краевых условий задачи. Например, книга [1] дает практически все о полиномах Чебышёва, что нужно при решении практических задач. Однако в примерах решения задач учет краевых условий там осуществляется путем подбора вручную уравнений для нужных коэффициентов. Похожий подход применялся в [3], p. 119.

Эта же проблема возникает в т.н. $\tau $-методе Ланцоша (см. [4]), который также является по сути спектральным. Учет краевых условий в $\tau $-методе проводится путем добавления уравнений для некоторых малых параметров ${{\tau }_{i}}$, $i = 1,2, \ldots $, а также выбрасывания некоторых уравнений для коэффициентов разложения решения из полученной переопределенной системы, но каждый раз в индивидуальном порядке.

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

Другой (и весьма малоизвестный) способ учета краевых условий был предложен Бабенко (см. [2]).

Бабенко предложил способ аппроксимации некоторых типовых дифференциальных операторов с учетом граничных условий так, как если бы аппроксимировались интегральные. Далее нужный дифференциальный оператор собирается из имеющихся прототипов так, что его граничные условия достаточно хорошо аппроксимируются (подробности см. в [5]).

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

Таким образом, проблема учета краевых условий в спектральных методах остается (с практической точки зрения) открытой и пока не имеет общего и эффективного решения.

В этой статье мы даем, в определенном смысле, полное решение этой задачи для конкретного спектрального метода – разложения функций на интервале $[0,1]$ по смещенным полиномам Лежандра (т.е. с единичным весом).

Выбор этих полиномов обусловлен тем, что разложения по ним также сходятся равномерно (при определенных условиях), т.е. полученные аппроксимации будут давать методы без насыщения.

Вторая причина, по которой мы выбрали эти полиномы – это наличие приложений данного метода к аналитической теории чисел, что возможно при выборе рациональной арифметики в расчетах и использовании компьютерной алгебры (CAS).

Учет краевых условий в разложениях по полиномам Чебышёва, на самом деле, осуществляется таким же образом. Однако эти разложения больше подходят для численных расчетов, и им будет посвящена отдельная статья.

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

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

2. АРИФМЕТИКА РАЗЛОЖЕНИЙ ПО ПОЛИНОМАМ ЛЕЖАНДРА

Везде в этой статье рассматриваются аналитические функции на отрезке $[0,1]$, которые могут иметь особенности на концах интервала. Эти функции раскладываются в ряды по смещенным полиномам Лежандра, и эти разложения сходятся (вообще говоря) равномерно на любом подынтервале отрезка $[0,1]$.

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

Пусть $\tilde {P}(n,x)$ – это обычные полиномы Лежандра, ортогональные на отрезке $[ - 1,1]$ с единичным весом. Мы используем смещенные полиномы Лежандра, $P(n,x) = \tilde {P}(n,2{\kern 1pt} x - 1)$, ортогональные на отрезке $[0,1]$,

Пусть $N$ – это размерность аппроксимации, т.е. любая функция $f(x)$ в этом разделе понимается как отрезок ее разложения по полиномам $P(n,x)$ до степени $N - 1$, т.е. $f \in {{\mathcal{P}}_{N}}$. Это представление функций аналогично представлению чисел в арифметике с фиксированной разрядной сеткой. Ясно, что дифференцирование и интегрирование таких функций нарушают это представление и требуют переразложения результата операции по полиномам $P(n,x)$ до номера $N - 1$ включительно. Это же справедливо для умножения функций. Например,

где равенство понимается в проективном смысле, т.е. при известных коэффициентах ${{a}_{n}}$ коэффициенты ${{b}_{n}}$ необходимо определить так, чтобы разница этих функций была ортогональна полиномам $P(n,x)$, $n < N$.

Таким образом, любая функция $f(x)$ представлена $N$-мерным вектором $f = {{\left\langle {{{a}_{0}},{{a}_{1}}, \ldots ,{{a}_{{N - 1}}}} \right\rangle }^{t}}$, а операции умножения на $x$, дифференцирование и т.д. – это линейные операторы в $N$-мерном векторном пространстве.

Мы будем обозначать векторы обычными, а линейные операторы – большими латинскими буквами. Введем следующие обозначения:

$e$ – единичный вектор, т.е. $e(x) = 1$, $e = {{\left\langle {1,0, \ldots ,0} \right\rangle }^{t}}$;

$E$ – единичная матрица;

$X$ – оператор умножения на $x$;

$X2$ – оператор умножения на ${{x}^{2}}$;

$x$ – вектор, соответствующий функции $x$, т.е. $x = {{\left\langle {\frac{1}{2},\frac{1}{2},0, \ldots ,0} \right\rangle }^{t}}$;

$D$ – оператор дифференцирования по $x$;

$I$ – оператор формального интегрирования по $x$.

Таким образом,

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

Заметим также, что $X.X \ne X2$, хотя эти две матрицы отличаются только в одном элементе, $(X2 - X.X)[N,N] = {{N}^{2}}{\text{/}}(4{\kern 1pt} {{N}^{2}} - 1){\text{/}}4$. Однако оператор $X2$ все же лучше в численных (но не символьных) расчетах, чем $X.X$, но не настолько, чтобы вводить отдельные операторы $X3,\;X4$ и т.д.

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

Матрицы $X$ и $X2$ являются соответственно трех- и пятидиагональными. Приведем для примера эти матрицы для размерности $N = 6$:

Для матрицы $X$ диагонали образуют последовательности, начиная с нижней, соответственно, Для матрицы $X2$ диагонали образуют последовательности, начиная с нижней, соответственно,

Проще всего устроена матрица дифференцирования $D$, а сложнее всего – матрица интегрирования $I$. Эти матрицы для размерности $N = 6$ имеют вид:

Элемент с индексами $[m,n]$ в верхнетреугольной матрице $D$ имеет вид Матрица $I$ является трехдиагональной, за исключением первой строки. Элементы главной диагонали нулевые. Ненулевой элемент с индексами $[m,n]$ в матрице $I$ имеет вид где функция номера столбца $c(n)$ дается формулой

Матрица $X$ размерности $N$ обладает важным свойством, которое нам понадобится. Справедливо

Предложение 1. Матрица $P(N,X)$ нулевая, $N = 1,2, \ldots $.

Доказательство. Введем вектор $p = {{\left\langle {P(0,x),P(1,x), \ldots ,P(N - 1,x)} \right\rangle }^{t}}$. Для трехдиагональной матрицы Якоби $J$ размерности $N$, ассоциированной с системой ортогональных полиномов Лежандра $\{ P(n,x),n = 0,1, \ldots \} $, справедливо тождество

которое является прямым следствием трехзвенного линейного рекуррентного соотношения для этих полиномов (см. [6]).

Результат умножения любого вектора $f \in {{\mathcal{P}}_{N}}$ на $x$ можно записать как скалярное произведение этого вектора с левой частью (1). Поэтому транспонируя (1), получим $X = {{J}^{t}}$.

Но спектр матрицы $J$ – это корни полинома $P(N,x)$, что следует из (1). Поэтому матрица $J$, а значит, и $X$ аннулируются своим характеристическим полиномом $P(N,x)$ (теорема Кэли–Гамильтона). Что и требовалось доказать.

Как следствие, получаем,

Также при нечетном $N$ все матрицы $P(M,X)$ для нечетных $M$ вырождены (но ненулевые). Это связано с тем, что у всех полиномов $P(M,x)$ нечетной степени есть общий рациональный корень $x = 1{\text{/}}2$. Для четных $N$ все матрицы $P(M,X)$, $M \ne N$ невырожденны.

Обратим внимание, что последняя строка матрицы дифференцирования $D$ всегда нулевая, т.е. этот оператор всегда вырожден, а его матрица имеет ранг $N - 1$. Это прямое следствие того, что ОДУ $y{\kern 1pt} '(x) = f(x)$ имеет одномерное многообразие решений. Иными словами, для однозначного определения решения $y(x)$ необходимо задать начальное (либо краевое) условие.

Эти условия можно учесть совершенно естественным образом в последней строке матрицы D. Поскольку

то последнюю строку матрицы $D$ надо заполнить этими значениями для учета, соответственно, начального (т.е. при $x = 0$) или краевого (т.е. при $x = 1$) условия (или их линейной комбинации). Затем нужно заменить последний элемент вектора $f$ в правой части на нужное начальное (краевое) условие. Тогда полученная линейная система будет невырожденна, а полученное решение $y(x)$ будет автоматически иметь нужное начальное (краевое) условие.

Последнее утверждение можно сформулировать точно. Справедливо

Предложение 2. Пусть правая часть $f(x)$ ОДУ $y{\kern 1pt} '(x) = f(x)$ имеет вид

т.е. последний элемент вектора $f = f(x)$ $($который заменяется на граничное условие$)$ равен нулю. Тогда решение $y(x)$ находится точно в виде полинома описанным выше способом при любом начальном $($краевом$)$ условии или их линейной комбинации.

Доказательство. Для того, чтобы ОДУ $y{\kern 1pt} '(x) = f(x)$ имело решение $y \in {{\mathcal{P}}_{N}}$ необходимо и достаточно, чтобы правая часть $f$ имела вид (2). Любое линейное краевое условие учитывается подбором константы интегрирования, для чего достаточно одного линейного уравнения. Что и требовалось доказать.

Можно дать еще одну интерпретацию описанного способа решения ОДУ $y{\kern 1pt} '(x) = f(x)$ с полиномиальной правой частью.

Как было отмечено, ранг матрицы аппроксимации оператора $D$ равен $N - 1$. Поэтому если решать полученную вырожденную линейную систему методом Гаусса с выбором главного элемента по всей матрице, то на последнем шаге $N$ получится уравнение $0 \times \tilde {y}[N] = \tilde {f}[N]$, где тильда символизирует операции, проделанные в методе Гаусса. Последнее уравнение может быть совместно, только если $\tilde {f}[N] = 0$, т.е. $f(x)$ – это полином вида (2).

Практически можно взять любое значение $\tilde {y}[N]$, а затем скорректировать начальное (краевое) условие полученного решения, как требуется. Эта операция полностью эквивалентна неопределенному интегрированию функции, т.е. третий способ решения ОДУ $y{\kern 1pt} '(x) = f(x)$ – это вычисление решения $y = I.f$, где вектор $f$ вида (2), с последующей корректировкой константы интегрирования.

В операторе ${{D}^{2}} = D.D$ будут, очевидно, две нулевые последние строки, которые следует заполнить нужными краевыми условиями. Вообще, для учета краевых условий нашим методом необходимо и достаточно знать следующие векторы-строки, линейные комбинации которых заполняют последние строки нужных линейных дифференциальных операторов. Эти векторы размерности $N$ имеют вид

где $m$ обозначает номер производной полинома $P(n,x)$ по $x$. Векторы ${{i}_{m}}$ отвечают за граничные значения функции при $x = 0$, а ${{b}_{m}}$ – за граничные значения при $x = 1$. Эти векторы устроены весьма просто: причем этот список может быть легко продолжен (с использованием CAS).

Любые начальные либо краевые линейные условия для функций $f(x)$, $f{\kern 1pt} '(x)$, $f{\kern 1pt} '{\kern 1pt} '(x)$ и т.д. и их комбинаций получаются в виде линейных комбинаций перечисленных векторов.

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

Приведем численный пример, на котором видны преимущества данного метода. Рассмотрим спектральную задачу на интервале $x \in [0,1]$:

с нестандартными краевыми условиями. Эта задача приводилась в [7] для демонстрации приложения $\tau $-метода.

Спектральная задача (3) решается в явном виде. Несложные вычисления показывают, что все собственные значения являются решениями трансцендентного уравнения $x + \tan x = 0$, т.е.

Если решать эту задачу методом Бабенко, то необходимо сперва сконструировать интерполяционный полином Чебышёва с указанными краевыми условиями. Затем этот полином нужно дважды продифференцировать и получить таким образом конечномерную спектральную задачу. Мы проделали эти вычисления, но они весьма громоздки, и мы их опускаем.

Согласно нашему методу, нужно взять матрицу $A = D.D$ и заменить ее две последние стороки на векторы ${{i}_{0}}$ и ${{b}_{0}} + {{b}_{1}}$ (в любом порядке). Полученная спектральная задача размерности $N = 10$ дает приближенные значения ${{a}_{1}},{\kern 1pt} \;{{a}_{2}},\;{\kern 1pt} {{a}_{3}}$, соответственно, с шестью, тремя и одним десятичными разрядами при вычислениях в обычной арифметике, что соответствует точности, полученной для этой задачи в [7] с помощью $\tau $-метода той же размерности.

При увеличении $N$ количество полученных собственных значений увеличивается, а их точность возрастает, т.е. спектр стабилизируется. Никаких “паразитных” собственных значений в нашем методе не появляется, т.е. данный метод является методом без насыщения.

В модификации $\tau $-метода, которая использовалась в [7], для учета краевых условий к разложению функции добавляются дополнительные члены, уравнения для коэффициентов которых получаются из краевых условий. Это требует весьма громоздких вычислений, предназначенных только для данного конкретного случая. В результате для системы (3) $\tau $-метод в [7] дает спектральную задачу для пучка матриц, т.е. $Ay = \lambda {\kern 1pt} By$, где матрицы $A$ и $B$ имеют блочную структуру, собранную из отдельных кусков.

Отметим также, что, помимо учета краевых условий, наш подход схож с операторным подходом к $\tau $-методу, предложенному в [8]. Однако в [8] используются бесконечные матрицы, которые обрезаются нужным образом для получения конечномерных аппроксимаций операторов. Затем эти конечномерные матрицы обрамляются дополнительными блоками (примерно как в [7]) для учета краевых условий.

3. КВАДРАТУРЫ ГАУССА И ПАДЕ-АППРОКСИМАЦИИ

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

Вычисления с расширенной разрядной сеткой, а также в рациональной арифметике, проделанные нами в CAS Maple, позволили установить связь нашего метода с квадратурами Гаусса и паде-аппроксимациями функций.

Напомним, что в классической $M$-точечной квадратуре Гаусса с единичным весом на интервале $[ - 1,1]$ требуется вычислить $M$ корней полинома Лежандра $\tilde {P}(M,x)$ и $M$ весов квадратурной формулы (см. [9]). На самом деле достаточно вычислить половину из каждого набора ввиду их симметрии, но проще вычислить их все сразу. Это возможно сделать только в плавающей арифметике. Тогда полученная квадратура будет точна на полиномах степени не выше $2{\kern 1pt} M - 1$.

Поскольку при нашем подходе функция приближается полиномами Лежандра степени не выше $N - 1$, то можно было ожидать, что данный метод будет соответствовать квадратуре Гаусса порядка $[N{\text{/}}2]$. Приведем пример, показывающий, что точность на самом деле значительно выше.

Рассмотрим задачу Коши

решением которой является функция $y(x) = \arctan x$.

Выберем для этого примера размерность аппроксимации $N = 10$.

Правая часть уравнения (4) может быть сразу записана в виде вектора $r = (E + X{{2)}^{{ - 1}}}.e$ (см. теорему ниже), т.е. мы получаем приближенное разложение рациональной функции $r(x)$ по смещенным полиномам Лежандра, просто обратив нужную матрицу. Обозначим точные коэффициенты разложения через ${{\tilde {r}}_{k}}$:

Тогда Поскольку то мы вычислили интеграл от правой части (4), даже не решая дифференциального уравнения. Имеем

Уместно сравнить полученную точность с той, которую можно достичь с помощью квадратур Гаусса при вычислении интеграла (6). Обозначим соответствующие квадратуры ${{G}_{M}}$. Тогда

т.е. наш метод соответствует (в данном примере) 11-точечной квадратуре Гаусса.

Теперь вычислим решение задачи Коши (4) согласно нашему методу. Возьмем матрицу $A = D$ и заменим последнюю (нулевую) строку на вектор ${{i}_{0}}$. Далее заменим последний (в данном случае ненулевой) элемент вектора $r$ на ноль, согласно начальному значению решения. Получаем конечномерную аппроксимацию задачи (4), $A.y = r$, откуда находим вектор $y = {{A}^{{ - 1}}}.r$. Таким образом, получаем приближенное разложение арктангенса

Подставив $x = 1$ в эту формулу, получаем рациональное приближение числа $\pi {\text{/}}4$, которое оказывается в точности равным значению ${{r}_{1}}$ в (7), вычисленному ранее.

Вычислим паде-аппроксимации арктангенса различной размерности в рациональной арифметике, подставим туда $x = 1$ и сравним c приближением числа $\pi {\text{/}}4$, полученным нами. Оказалось, что

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

Коэффициенты ${{y}_{k}}$ отличаются от точных значений разложения арктангенса примерно так же, как в (5). На фиг. 1 показаны графики разности $\arctan x$ с их аппроксимациями полиномами Лежандра для $N = 10$. На фиг. 1 слева использовались точные (иррациональные) коэффициенты разложения, а на фиг. 1 справа – полученные нами.

Обратим внимание, что глобальная (или равномерная) погрешность нашей аппроксимации примерно соответствует ожидаемой для данной размерности (и даже лучше). Однако краевое значение в задаче (4) находится с точностью, превышающей равномерную на несколько порядков.

Последнее обстоятельство также не является случайным. Чтобы это увидеть, решим задачу (4) нашим методом еще раз, но без использования матрицы $X2$. Сохраним те же обозначения, но возьмем правую часть уравнения (4) в виде вектора $r = (E + X.X{{)}^{{ - 1}}}.e$. Тогда

т.е. точность аппроксимации правой части (4) несколько хуже, чем ранее. Это же справедливо для нового вектора $y$. Однако теперь что соответствует 10-точечной квадратуре Гаусса (см. выше).

Данное соответствие нашего рационального приближения интеграла его квадратуре Гаусса на самом деле является полным совпадением.

Возьмем произвольную рациональную функцию $r(x)$ с рациональными коэффициентами (можно и комплексно-рациональными или даже символьными). Рассмотрим задачу Коши

Ее решением является функция $z(x) = \int_0^x \,r(t){\kern 1pt} dt$, если она существует. Но нас интересуют формальные символьные вычисления в конечномерной аппроксимации.

Для задачи (9) в случае общего положения (см. ниже) вычислим следующие массивы для размерностей аппроксимации $N = 1,2, \ldots $

${{S}_{1}}[N] = {{r}_{1}}$, где вектор $r = {{\left\langle {{{r}_{1}}, \ldots ,{{r}_{N}}} \right\rangle }^{t}}$ вычисляется по формуле

где коммутирующие матрицы $p(X)$ и $q(X)$ получаются подстановкой матрицы $X$ вместо независимой переменной $x$ в полиномы $p(x)$ и $q(x)$. Напомним, что $X2$ мы теперь не используем.

${{S}_{2}}[N] = \sum\nolimits_{k = 1}^N {{{z}_{k}}} $, где вектор $z$ вычисляется так же, как вектор $y$ для (4). То есть берем матрицу $A = D$, заменяем последнюю (нулевую) строку на вектор ${{i}_{0}}$. Далее заменяем последний (ненулевой) элемент вектора $r$ на ноль, т.е. берем $r[N] = 0$, согласно начальному значению решения. Затем вычисляем вектор $z = {{A}^{{ - 1}}}.r$, который дает приближенное разложение функции $z(x)$, если она существует. Иначе это просто формальный вектор. Отметим, что ${{S}_{2}}[1] = 0$ в силу способа вычисления.

Этот массив можно также вычислять альтернативным способом с помощью оператора формального интегрирования, как это описано в разд. 2.

${{S}_{3}}[N] = {{G}_{N}}(r)$, т.е. ${{S}_{3}}[N]$ – это значение N-точечной квадратуры Гаусса функции $r(x)$ на интервале $[0,1]$.

Напомним, что корни ${{x}_{n}}$ полиномов Лежандра $P(N,x)$ могут быть вычислены (в принципе) в виде радикалов только до $N = 9$ включительно. Это же относится к весам ${{w}_{n}}$ квадратурной формулы Гаусса,

Разумеется, практически делать это бессмысленно, поэтому массив ${{S}_{3}}[N]$ (в отличие от предыдущих) дается в плавающей арифметике, но с любым (в принципе) числом десятичных разрядов.

Наложим очевидное ограничение на полюса функции $r(x)$, т.е. если $q(x)$ делится на полином Лежандра $P(M,x)$, то в массивах ${{S}_{1}},\;{{S}_{2}},\;{{S}_{3}}$ следует пропустить размерность $N = M$ (см. предложение 1). А в случае если $q(x)$ имеет корень $x = 1{\text{/}}2$, то эти массивы вычисляются только для четных $N$.

Заметим, что функция $r(x)$ не может иметь полюс ($ \ne \,1{\text{/}}2$) в одном из корней полинома $P(M,x)$ и не иметь полюсов во всех остальных корнях.

С учетом этих условий справедлива

Теорема. Массивы ${{S}_{k}}$, $k = 1,2,3$, совпадают.

Доказательство. Напомним, что $x = X.e$ по определению, т.е. вектор $x \in {{\mathcal{P}}_{N}}$, соответствующий функции $x$, получается применением матрицы $X$ к вектору $e$. Аналогично получается вектор $p = p(X).e$, соответствующий любому полиному $p(x)$. Но эти же рассуждения применимы для общей функции от матрицы, $f(X)$, т.е. вектор $f = f(X).e$ – это коэффициенты приближенного разложения функции $f(x)$ по смещенным полиномам Лежандра при условии, что матрица $f(X)$ определена.

Функция от матрицы $f(X)$ однозначно определяется значениями функции $f(x)$ на спектре матрицы $X$ (см. [10]), т.е. значениями этой функции в корнях полинома $P(N,x)$ (так как матрица $X$ диагонализируема, т.е. производные $f(x)$ не нужны).

Для рациональной функции $r(x)$ не требуется вычислять ее значения на спектре матрицы $X$, а можно формально подставить матрицу $X$ вместо переменной $x$ (см. [10]). Поэтому совпадение массивов ${{S}_{1}}$ и ${{S}_{2}}$ следует из совпадения двух интегралов.

Величина ${{S}_{1}}[N]$ – это первая компонента вектора правой части (9), т.е. нулевой коэффициент ее аппроксимации полиномами $P(n,x)$. Иными словами, это аппроксимация интеграла $r(x)$ на отрезке $[0,1]$. Величина ${{S}_{2}}[N]$ получается из аппроксимации интеграла $\int_0^x \,r(t){\kern 1pt} dt$ как функции на отрезке $[0,1]$ подстановкой $x = 1$ в эту функцию, т.е. ${{S}_{1}} = {{S}_{2}}$.

Вторая часть утверждения теоремы, т.е. ${{S}_{1}} = {{S}_{3}}$, следует из представления матричной функции $r(X)$ интерполяционным полиномом Лагранжа–Сильвестра (см. [10]). Имеем

где $\{ {{x}_{n}},{\kern 1pt} n = 1, \ldots ,N\} $ – это спектр матрицы $X$ (корни полинома $P(N,x)$), а ${{l}_{{N,n}}}(x)$ – это фундаментальные полиномы лагранжевой интерполяции.

Поскольку ${{l}_{{N,n}}}(x) \in {{\mathcal{P}}_{N}}$, то

где ${{w}_{n}}$ – это веса (10) квадратурной формулы Лежандра–Гаусса (см. [9]). Таким образом, вектор правой части (9) представим в виде поэтому ${{S}_{1}} = r(X).e[1]$ – это квадратура Лежандра–Гаусса функции $r(x)$. Что и требовалось доказать.

Следствие. Любая квадратура Лежандра–Гаусса на интервале $[0,1]$ от любой рациональной функции $r(x)$ $($с учетом сделанных оговорок$)$ является рациональным числом, т.е.

Этот результат может быть (относительно легко) проверен в CAS до размерности $N = 5$ и, в принципе, до размерности $N = 9$, хотя выражение для корня полинома $\tilde {P}(8,x)$ занимает страницу. Для $N > 9$ этот рузультат может быть проверен только численно, но зато с любой, в принципе, точностью.

Возьмем, например, $r(x) = 1{\text{/}}x$, для которой решение задачи Коши (9) $z(x)$ не существует. Но это не мешает нам вычислить квадратуры Гаусса, так как все корни ${{x}_{n}}$ лежат внутри интервала. Получим в точности числа $2{\kern 1pt} H(N)$, $N = 1,2, \ldots $, где $H(N)$ – это гармоническое число.

Возьмем теперь $r(x) = q{\text{/}}(1 + q{\kern 1pt} x)$. Тогда $z(x) = \ln (1 + q{\kern 1pt} x)$ является формальным решением задачи Коши (9). По нашему алгоритму получим для $N = 1,2, \ldots $,

Квадратуры Гаусса (11) получаются весьма громоздкими и в плавающей арифметике, но при подстановке численных значений параметра $q$ совпадают с рациональными функциями (12).

Формула (12) также показывает, почему полюс $x = 1{\text{/}}2$ при $q = - 2$ неприемлем. Он соответствует корню $x = 0$ полиномов Лежандра нечетной степени. Других рациональных корней у полиномов Лежандра не существует.

Предложение 3. Последовательность рациональных приближений $(12)$, т.е. ${{S}_{1}}[N]$, $N = 1,2, \ldots $, совпадает с диагональными паде-аппроксимациями ${{\left. {P[N,N](z(x))} \right|}_{{x = 1}}} = P[N,N](\ln (1 + q))$.

Доказательство. Докажем другое утверждение, которое равносильно этому согласно доказанной теореме. Последовательность квадратур Гаусса функции $r(x) = q{\text{/}}(1 + q{\kern 1pt} x)$, т.е. ${{S}_{3}}[N]$, $N = 1,2, \ldots $, совпадает с этими паде-аппроксимациями. В самом деле,

в силу линейности квадратуры. Но квадратуры Гаусса точны на полиномах степени $M \leqslant 2{\kern 1pt} N - 1$, поэтому Слева здесь стоит рациональная функция от $q$ нужного вида, т.е. отношение полиномов степени $N$. Что и требовалось доказать.

Например, подставив $q = 1$ в (12), получаем последовательность рациональных приближений для $\ln 2$:

которая совпадает с последовательностью, полученной в [11] с помощью $\varepsilon $-алгоритма Винна (Wynn). Пользуясь случаем, исправляем там опечатку: ($8{\text{/}}13 \to 9{\text{/}}13$).

Эта последовательность соответствовала в [11] поддиагональной паде-аппроксимации, так как там рассматривалась функция $\ln (1 + x){\text{/}}x$. Отметим также, что в [11] было выбрано неудачное преобразование эквивалентности, в результате чего полученная цепная дробь выглядела сложнее, чем она есть (см. разд. 4).

Совпадение паде-аппроксимаций функции, представленной рядом Стилтьеса, т.е. где коэффициенты разложения – это моменты интеграла Стилтьеса, с квадратурами Гаусса не является новым результатом (см. [12]) (авторы там перепутали $[n - 1,n]$- и $[n,n - 1]$-аппроксимации).

Таким образом, существенно нелинейная операция – вычисление диагональной паде-аппроксимации $\ln (1 + x)$ с помощью $\varepsilon $-алгоритма Винна сводится нашим методом к линейным операциям с матрицами.

Также теперь ясно, почему не удалось получить точную паде-аппроксимацию для арктангенса. Дело в том, что наш алгоритм линеен, поэтому мы получили сумму двух различных паде-аппроксимаций. А именно,

т.е. последовательность $\{ {{S}_{1}}[N]\} $ соответствует сумме двух паде-аппроксимаций,

Хотя любая рациональная функция представима в виде суммы простых дробей, следствие, которое мы вывели из теоремы, не следует из совпадения квадратур Гаусса с паде-аппроксимациями простых дробей по очевидным причинам.

Как показано в [11], из последовательности (13) следует иррациональность числа $\ln 2$ по критерию Дирихле, т.е.

где $p{\text{/}}q$ – это рациональные приближения в (13). Этим свойством обладают некоторые (но далеко не все) рациональные последовательности, полученные с помощью паде-аппроксимаций, и мы займемся этим в следующем разделе.

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

Оказалось, что если эти коэффициенты исходно рациональны, то наш метод в ряде случаев вычисляет их точно в рациональной арифметике. Рассмотрим уравнение

решением которого является функция $y(x) = \ln x$ при выборе краевого условия $y(1) = 0$. Это также первый пример, где мы рассматриваем более сложный оператор, чем просто дифференцирование. Однако способ решения этого уравнения точно такой же.

Берем матрицу $A = X.D$ и заменяем последнюю строку на вектор ${{b}_{0}}$. Последний элемент вектора $r = e$ здесь ноль. Вычисляем вектор $y = {{A}^{{ - 1}}}.r$, который дает приближенное разложение $\ln x$. Однако в данном случае все коэффициенты разложения, кроме последнего, $y[N]$, оказываются точными. Таким образом, получаем

Если бы мы решали уравнение $y{\kern 1pt} '(x) = 1{\text{/}}x$ этим способом, что возможно, так как матрица $X$ обратима, то мы получили бы почти такой же результат, т.е. все коэффициенты находятся правильно, кроме (по неизвестным пока причинам) предпоследнего.

Приведем еще один пример. Рассмотрим уравнение

общим решением которого является функция $y(x) = x{\kern 1pt} (C + \ln x)$. Поэтому возьмем краевое условие $y(1) = 0$ и решим эту задачу как предыдущую. Тогда при каждом $N$ мы получим все коэффициенты разложения $y(x) = x{\kern 1pt} {\kern 1pt} \ln x$, кроме последнего, $y[N]$. Таким образом, получаем

Полученные разложения могут быть использованы для решения неалгебраических ОДУ нашим методом (см. разд. 5). Хотя равномерная аппроксимация решений окажется при этом соответствующей ожидаемой, т.е. (вообще говоря) не может быть лучше, чем аппроксимация $\ln x$ полиномами Лежандра, но аппроксимация решений в концах интервала может быть очень хорошей, что мы уже наблюдали (см. фиг. 1).

4. ЦЕПНЫЕ ДРОБИ И ИРРАЦИОНАЛЬНОСТЬ

Как известно, быстро сходящиеся последовательности рациональных чисел играют ключевую роль в доказательствах иррациональности их предела (см. [11] и ссылки там).

При этом, как оказалось во всех без исключения случаях (а их можно пересчитать по пальцам одной руки), если последовательность $\{ {{p}_{n}}{\text{/}}{{q}_{n}}\} $ дает иррациональность своего предела $x$ (например, константы Апери) по критерию Дирихле, то она сходится к $x$ быстрее (в обычном, численном смысле), чем обыкновенная цепная дробь числа $x$. А значит (см. [11]), эта последняя может быть однозначно восстановлена по последовательности $\{ {{p}_{n}}{\text{/}}{{q}_{n}}\} $ с помощью функций $K(x,n)$, которые мы определили в [11].

Иными словами, в этом случае последовательность $K({{p}_{n}}{\text{/}}{{q}_{n}},n)$ совпадает с последовательностью подходящих дробей числа $x$, т.е. $K({{p}_{n}}{\text{/}}{{q}_{n}},n) = K(x,n)$, $n > {\text{const}}$.

Но это последнее свойство может также выполняться, если критерий Дирихле для последовательности $\{ {{p}_{n}}{\text{/}}{{q}_{n}}\} $ не выполняется (см. [11]). Например, для последовательности $\{ {{S}_{1}}[N] = {{r}_{1}}\} $ рациональных приближений числа $\pi {\text{/}}4$, которую мы получили в задаче (4), $K({{S}_{1}}[N],N) = K(\pi {\text{/}}4,N)$, $N = 1,2, \ldots $, т.е. мы получаем все подходящие дроби числа $\pi {\text{/}}4$. При этом, например,

Однако наибольший интерес все же представляют последовательности, в структуре которых удается выявить некоторые закономерности. Этим мы и займемся в этом разделе.

Последовательности диагональных, а также над- или поддиагональных паде-аппроксимаций и их комбинации ($qd$-алгоритм) некоторых функций могут давать очень хорошие рациональные приближения иррациональных чисел. Например, последовательность (13) дает иррациональность числа $\ln 2$ по критерию Дирихле (см. [11]).

Но у нас также имеется более общая последовательность (12). Как оказалось, эта последняя значительно более удобна для выявления закономерностей в ее структуре, чем ее частный вариант при $q = 1$. Это связано с тем, что параметр $q$ не дает CAS приводить дроби к наименьшему знаменателю, что в данном случае контрпродуктивно.

Закономерность в последовательности (12) становится видна, если преобразовать ее в цепную дробь Эйлера (по формулам Д. Бернулли). Обозначим (12) (заменив там $q \to t$) как последовательность $s(n) = {{p}_{n}}{\text{/}}{{q}_{n}}$ подходящих дробей, где ${{p}_{n}}$ и ${{q}_{n}}$ – это, соответственно, числитель и знаменатель дробей (12). Вычислим

где ${{p}_{0}} = 0$ и ${{q}_{0}} = 1$. Тогда или где закономерность в частных дробях очевидна.

Хотя нам не удалось получить Паде-аппроксимации арктангенса с помощью нашего алгоритма напрямую, но применение формулы

в (16) дает такое разложение. После некоторых преобразований получаем известное разложение

Применим наш подход к некоторым ОДУ, зависящим от параметра. Рассмотрим задачу Коши

где $t$ – это независимый параметр. Это также пример, как формально однородное уравнение решается как неоднородное за счет начального условия, т.е. вектор в правой части $r = 0$, но затем $r[N] = 1$.

Применив наш алгоритм и преобразование (15), получим

т.е. после начального куска цепная дробь становится “квазипериодической”.

Однако в данном случае мы получили феномен, который ранее (насколько нам известно) не наблюдался. То есть мы получили разложение, которое (при $t = 1$) содержит обыкновенную цепную дробь некоторой константы. Более того, дробь (17) (при $t = 1$) сходится на несколько порядков быстрее, чем обыкновенная цепная дробь константы ${\mathbf{e}}$.

Данный феномен объясняется довольно просто. Избавившись от начального куска дроби (17) и достроив ее сверху согласно очевидным закономерностям, получим известное разложение

Таким образом, преобразование Мёбиуса с рациональными коэффициентами иррационального числа может радикальным образом изменить скорость сходимости обыкновенной цепной дроби новой константы по сравнению с исходной. Например,

Это свойство могло бы найти применение в доказательствах иррациональности, но (насколько нам известно) никогда не применялось.

Применим теперь формулу

в (18). Тогда (после преобразования эквивалентности) получим классическое разложение для тангенса:

Это разложение Ламберт и Лежандр использовали в доказательстве иррациональности чис-ла π (см. [13], p. 289–296). А именно, согласно (очень простому) критерию иррациональности Ламберта–Лежандра формула (19) может давать только иррациональные числа для рациональных $x$. Но $\tan (\pi {\text{/}}4) = 1$, поэтому $\pi $ не может быть рациональным.

По тем же причинам из дроби (18) следует иррациональность $\exp (r)$ для $0 \ne r \in \mathbb{Q}$. Поэтому и $\ln (r)$ иррационально для $1 \ne r \in \mathbb{Q}$.

Заметим, что основную часть доказательства Лежандра в [13] занимал как раз вывод цепной дроби (19). Лежандр использовал для этого функции Бесселя, причем, по-видимому, до самого Бесселя. А мог бы, как теперь ясно, применить свои полиномы.

Приведем еще один пример. Рассмотрим задачу Коши

где $s$ и $t$ – это независимые параметры. Применив наш алгоритм и преобразование (15), получим т.е., как и раньше, после начального куска цепная дробь становится “квазипериодической”. Избавившись от начального куска дроби (20) и достроив ее сверху согласно очевидным закономерностям, получим разложение которое весьма напоминает (18).

Применим теперь тождество ${{i}^{{ - 2{\kern 1pt} i}}} = {{{\mathbf{e}}}^{\pi }}$, т.е. сделаем подстановку $s = - 2{\kern 1pt} i$, $t = i - 1$ в (21). Тогда получим (после преобразования эквивалентности)

Цепная дробь для ${{{\mathbf{e}}}^{\pi }}$ была получена недавно в [14]. Причем, как полагает автор в [14], это было сделано впервые. Однако наша дробь устроена проще.

Заметим также, что цепная дробь (22) сходится к своему пределу медленнее, чем его подходящие дроби, хотя число ${{{\mathbf{e}}}^{\pi }} = ( - {{1)}^{{ - i}}}$ трансцендентно (как алгебраическое число в степени мнимая квадратичная иррациональность (см. [15] и ссылки там)).

5. КОНСТАНТЫ ЭЙЛЕРА И ЭЙЛЕРА–ГОМПЕРТЦА

Как было отмечено в [16], константа Эйлера $\gamma $ и константа Эйлера–Гомпертца $\delta $ имеют весьма похожие интегральные представления. В этом разделе мы продолжим эту аналогию.

Рассмотрим функцию

где ${\text{Ei}}()$ – это интегральная экспонента. Функция $y(x)$ удовлетворяет ОДУ: общим решением которого является функция

Легко проверить, что $y(0) = - \gamma + C$ и $y(1) = y{\kern 1pt} '(1) = \delta + C\exp (1)$. Поэтому, к сожалению, из уравнения (23) нельзя получить сразу обе константы, $\gamma $ и $\delta $, а только их линейную комбинацию с константой ${\mathbf{e}}$. Иными словами, задача Коши для уравнения (23), $y(0) = 0$, дает рациональные приближения для константы

Эта формула следует из представления (см. [17]) где $F()$ – это гипергеометрическая функция, т.е., как легко проверить,

Задача Коши для уравнения (23) решается точно так же, как мы уже делали неоднократно. Отличие состоит в том, что эта задача сингулярна ($y{\kern 1pt} '(0) = + \infty $), а также что ОДУ не является алгебраическим.

Таким образом, решается линейная алгебраическая задача $A.y = r$, где $A = D - E$ с последующей корректировкой последней строки, а вектор $r$ получается с помощью формулы (14), и $r[N] = 0$, согласно начальному условию.

Мы проделали эти вычисления в рациональной арифметике до размерности аппроксимации $N = 100$. Полученная в результате последовательность рациональных приближений $\{ s(n),n = 1, \ldots ,100\} $ константы $\eta $ не дает иррациональности этой константы по критерию Дирихле, но приближает ее лучше в численном смысле, чем последовательность подходящих дробей константы $\eta $, т.е., согласно [11], $K(\eta ,n) = K(s(n),n)$, $n > 11$. Иными словами, получена последовательность подходящих дробей константы $\eta $ начиная с 12-й. Например,

Здесь мы воспользовались модификацией алгоритма, т.е. замедлили скорость сходимости полученной последовательности рациональных приближений с помощью функций $K(x,n)$, введенных в [11].

Данный модифицированный алгоритм (т.е. вычисление $n$-й подходящей дроби от $n$-го рационального приближения) обладает сертификатом качества (см. [11]). Он гарантированно дает подходящие дроби нужной константы либо не дает, о чем сигнализирует. При этом не надо знать численного значения данной константы, а все вычисления проводятся в рациональной арифметике.

В случае если алгоритм “дает сбой”, т.е. очередное рациональное приближение оказывается хуже, чем приближение подходящей дробью этой константы (что может быть, например, если в цепной дроби встретился нетипично большой знаменатель), то это никак не влияет на последующие подходящие дроби, если данный член просто опустить. Такой случай встречался в разложении $\zeta (6)$ (см. [11]). Полученная по формулам (15) цепная дробь просто окажется некоторым сжатием исходной.

Таким образом, свойство последовательности рациональных чисел $\{ s(n)\} $ сходиться к своему пределу $x$ быстрее (в численном смысле), чем последовательность подходящих дробей числа $x$, $\{ K(x,n)\} $, является внутренним свойством данной последовательности и проверяется в рациональной арифметике независимо от знания ее предела $x$.

Воспользуемся этими соображениями для проверки скорости сходимости ряда (24). Этот ряд весьма напоминает ряд Энгеля (см. [18]), а такие ряды обычно сходятся к своему пределу быстрее, чем подходящие дроби этого предела (см. [11]).

Гармоническое число в числителе (24) не должно сильно повлиять на скорость сходимости, так как оно растет как $\ln n$, что пренебрежимо мало по сравнению с факториалом.

Как оказалось, последовательность частичных сумм ряда (24), $\{ h(m)\} $, и в самом деле стремится к $\eta $ быстрее, чем подходящие дроби этой константы, т.е. $K(\eta ,m) = K(h(m),m)$, $m > 40$. Например,

В завершение этого примера решим задачу Коши $y(0) = - \gamma $ для уравнения (23), но на сей раз в символьном виде. То есть $\gamma $ и $\delta $ теперь просто символы. Найдем последовательность приближений для $y(1) = \delta $, $N = 1,2, \ldots $. Эта последовательность имеет вид $\tilde {\eta }(N) - {\mathbf{\tilde {e}}}(N){\kern 1pt} \gamma $, где тильда означает, что это рациональные приближения соответствующих констант. Они оказались в точности равны полученным ранее, т.е. в формулах (25) и (17).

Таким образом, если иметь последовательность хороших рациональных приближений для $\gamma $ (или $\delta $), то по данному алгоритму можно получить (почти) столь же хорошую последовательность рациональных приближений для $\delta $ (или $\gamma $).

6. СУММИРОВАНИЕ ГОЛОНОМНЫХ ПОСЛЕДОВАТЕЛЬНОСТЕЙ

Как мы видели, наш метод дает эффективные рациональные приближения коэффициентов разложения Фурье–Лежандра для функций, которые являются решениями линейных однородных ОДУ с полиномиальными коэффициентами. Такие функции называются голономными (см., например, [19]). Те неоднородные ОДУ, которые мы рассматривали, очевидно приводятся к однородным.

Оригинальный $\tau $-метод Ланцоша также предназначался для этого класса уравнений (см. [4]).

Последовательность чисел $\{ {{a}_{n}},n = 0,1, \ldots \} $ называется голономной, если она удовлетворяет линейному разностному уравнению с полиномиальными (от $n$) коэффициентами.

Если в качестве функции понимать формальный степенной ряд, то каждой голономной функции соответствует голономная последовательность чисел – коэффициентов ее формального разложения, а также наоборот: каждой голономной последовательности соответствует голономная функция – обыкновенная производящая функция этой последовательности (см. [19, p. 140]).

Хотя многие математики, продолжая традиции 19-го столетия, полагают, что “не каждый формальный ряд соответствует аналитической функции” (см. [19, p. 29]), однако это не так. Каждый формальный ряд соответствует аналитической функции, правда, возможно, не одной (теорема Бореля–Ритта, см. [20, p. 43]). Все зависит от интерпретации слова “соответствует”.

Возьмем ряд из [19, p. 29] в качестве примера. Мы использовали этот ряд в [21] как пример “голономного” суммирования последовательности. Рассмотрим формальный ряд

и найдем для него аналитическую функцию $y(x)$, такую, что $y(1) = 1$.

Формальный голономный ряд (26) удовлетворяет дифференциальному уравнению

которое неголономно, но приводится к нужному виду. Функция где $C$ – это константа интегрирования, является общим аналитическим решением этого уравнения.

Для любой константы $C$ эта аналитическая функция “соответствует” своему формальному разложению (26), которое является асимптотическим в полуплоскости ${\text{Re}}(x) > 0$, с трансфинитной экспоненциально малой добавкой.

Это означает, что $y(1)$, например, может принимать любое значение. Если взять

то получим $y(1) = 1$.

Чтобы получить разложение этой функции по полиномам Лежандра нашим методом, нужно взять матрицу $X2.D - E + X$ и заменить ее последнюю строку на вектор ${{b}_{0}}$ (см. разд. 2). Вектор правой части – это $r = - e$ с подстановкой $r[N] = 1$, согласно выбранному краевому значению. Заметим, что $X2.D = X.X.D$. Для $N = 10$ получим

что весьма близко к верхней равномерной оценке погрешности для этой функции на интервале $[0,1]$.

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

Таким образом, мы имеем способ суммирования голономных последовательностей и рядов, причем необязательно сходящихся в обычном смысле, который сводится к решению соответствующих голономных ОДУ и подстановке $x = 1$ в полученное решение.

Здесь наш метод приближенного решения голономных ОДУ может быть весьма полезен, так как эти уравнения, вообще говоря, неинтегрируемы. Но даже в случае их интегрируемости в квадратурах (с применением CAS) эти квадратуры могут быть очень сложно устроены и (почти) совершенно бесполезны.

Вычислим для примера константу Эйлера–Гомпертца путем суммирования расходящегося ряда (26) при $x = - 1$. Это классическая задача, которая решалась многими способами. Ланцош ее также использовал для демонстрации эффективности $\tau $-метода (см. [4]).

Сделаем замену $x \to - x$ в уравнении (27), затем в полученном уравнении сделаем замену $y(x) \to y(x){\text{/}}x$. Тогда получим сингулярную задачу Коши $y(0) = 0$ для уравнения

Мы выбрали такую форму записи уравнения для демонстрации одного эффекта, который встречается при решении сингулярных задач нашим методом. В данном примере конечномерная аппроксимация строится особенно просто: берется матрица $A = X2.D + E$ и правая часть $r = X.e$. Никакой корректировки матрицы $A$ и вектора $r$ делать не нужно, так как единственное ограниченное решение задачи (28) имеет начальное значение $y(0) = 0$ (т.е. $C = 0$ в (28)). В результате для размерностей аппроксимации $N = 10$ и $N = 20$ получим

Полученная последовательность рациональных приближений $\{ y(1)\} $ для $N = 1,2, \ldots $ сходится к $\delta $ весьма быстро, но все же медленнее, чем последовательность подходящих дробей $K(\delta ,N)$. Однако классическая цепная дробь

которая соответствует диагональной паде-аппроксимации ряда (28), ${{p}_{N}} = P[N,N](y(x)){{{\text{|}}}_{{x = 1}}}$, сходится к $\delta $ намного медленнее. Например,

Приведем два примера суммирования для очень медленно (и очень плохо) сходящихся голономных последовательностей. Рассмотрим ряды

где $L(n,x)$ – это полиномы Лагерра.

Второй из этих рядов мы вывели (без подробного доказательства) в [11] с помощью ряда Ньютона. Первый ряд в (29) получается из разложения $\ln x$ по полиномам Лагерра.

Предложение 4. Формулы (29) правильны.

Доказательство. Мы дадим “computer assisted proof”, что является стандартной практикой при работе с голономными рядами. Иными словами, средства CAS привлекаются для вывода формул, но сами формулы (в принципе) проверяемы вручную.

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

или Тогда можно проверить, что $u(x)$ и ${v}(x)$ удовлетворяют уравнениям с начальными условиями, определяемыми степенными разложениями $u(x)$ и ${v}(x)$. Эти уравнения интегрируемы. Имеем откуда следует $u(1) = \gamma $ и ${v}(1) = \delta $. Что и требовалось доказать.

Суммирование рядов (29) с помощью стандартных ускорителей сходимости (см., например, [22]), по-видимому, невозможно. Чтобы это увидеть, обозначим частичные суммы рядов (29) как $u(n)$ и ${v}(n)$ соответственно. Пользуясь известным рекуррентным соотношением для полиномов Лагерра, получаем $u(1) = 0$, $u(2) = 1{\text{/}}4$, $u(3) = 17{\text{/}}36$ и

а также ${v}(1) = 1$, ${v}(2) = 1$, ${v}(3) = 5{\text{/}}6$ и

На фиг. 2 показаны графики частичных сумм $u(n)$ и ${v}(n)$ до $n = 200$. Из фиг. 2 видно, что эти суммы осциллируют с медленно уменьшающейся амплитудой и увеличивающимся периодом. Такие ряды не могут суммироваться обычными средствами.

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

Мы уже неоднократно описывали алгоритм решения таких уравнений. Уравнения (31) решаются точно так же. В результате для размерностей аппроксимации $N = 10$ и $N = 20$ получим для $\gamma $:

а также для $\delta $:

Таким образом, хотя голономные последовательности для $\gamma $ и $\delta $ выглядят почти идентично (см. (29)), сам вид уравнений (31), а также результаты аппроксимации их решений (32) свидетельствуют, что константа $\gamma $ имеет более сложную природу в некотором “голономном” смысле. На это же указывает тот факт, что для рациональных приближений константы $\delta $ существуют линейные формы второго порядка, а для $\gamma $ пока что только третьего (см. [16]).

В заключение заметим, что наши методы обобщаются также на другие системы ортогональных полиномов, например, на полиномы Чебышёва. Последние гораздо удобнее для численного решения голономных ОДУ, для итерационного решения нелинейных уравнений, а также для решения разностных уравнений, которые возникают при обобщенном суммировании рядов (см. [21]). Однако это является темой отдельной работы.