Журнал вычислительной математики и математической физики, 2022, T. 62, № 2, стр. 199-216

1. ВВЕДЕНИЕ

Данная статья продолжает работу [1], где некоторые важные вопросы были оставлены для дальнейшего исследования.

В статье [1] рассматривалась задача интерполяции итерационных последовательностей, что можно интерпретировать как восстановление непрерывной динамической системы (НДС) по ее дискретной редукции. В самом деле, любую итерационную последовательность можно рассматривать как дискретную динамическую систему (ДДС), т.е. дискретное отображение (возможно, с запаздыванием). Таким образом, в [1] рассматривалась задача, в некотором смысле обратная той, которая возникает при редукции НДС к дискретным системам.

За последние примерно полвека интерес к ДДС значительно возрос во многом благодаря знаменитой статье Хенона и Хейлеса (см. [2]). Видимо, начиная с этой публикации ДДС стали рассматриваться как весьма содержательные модельные задачи, которые выражают самое существенное в динамике НДС.

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

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

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

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

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

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

Наряду с этим существует ряд работ (см. [4] и ссылки там), посвященных доказательству существования аналитических инвариантных кривых в некоторых модельных ДДС, как консервативных, так и диссипативных.

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

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

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

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

где F(y) является голоморфной функцией в некоторой окрестности нуля.

Отображение (1) дает пример простейшей ДДС ${{y}_{{n + 1}}} = F({{y}_{n}}) = {{F}^{{(n + 1)}}}({{y}_{0}})$, $n \in \mathbb{Z}$, $y \in \mathbb{R}$ или $y \in \mathbb{C}$. Как известно, формальные ряды непрерывной итерации (НИ) ${{F}^{{(x)}}}(y)$ (см. разд. 2) почти всегда расходятся, если $x \notin \mathbb{Z}$ (см. [1] и ссылки там).

Однако численные эксперименты в СХХ, а также со стандартным отображением Чирикова и в других задачах (см. [5]) показывают, что приближенные инвариантные кривые “идентичны натуральным”, т.е. ведут себя так же, как будто они являются аналитическими. Этот факт неоднократно упоминается в литературе в том плане, что численно установить отличие инвариантных кривых от аналитических не представляется возможным (см. [2]).

Это же наблюдение можно сделать при изучении НИ ${{F}^{{(x)}}}(y)$. Как отмечалось в [1], НИ ${{F}^{{(x)}}}(y)$, $x \in \mathbb{R}$, хорошо интерполируют дискретную последовательность ${{y}_{{n + 1}}} = F({{y}_{n}})$, $n \in \mathbb{Z}$, несмотря на заведомую расходимость рядов ${{F}^{{(x)}}}(y)$, $x \notin \mathbb{Z}$.

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

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

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

2. ИНВАРИАНТНЫЕ КРИВЫЕ НЕПРЕРЫВНОЙ ИТЕРАЦИИ

Напомним некоторые факты о НИ отображения (1), которые нам понадобятся (см. [1] и ссылки там).

Рассмотрим обычную композицию (итерацию) отображения (1)

Очевидно, что

НИ ${{F}^{{(x)}}}(y)$ обобщает это полугрупповое свойство на произвольные комплексные числа $x \in \mathbb{C}$, т.е. порождает непрерывную группу отображений.

НИ отображения (1) можно записать в виде

где ${{a}_{n}}(x)$ – однозначно определяемые полиномы степени не выше n.

Другая форма записи НИ использует генератор ${{g}_{1}}(y) = \omega (y)$ указанной группы отображений (см. [6]), т.е.

Вычисление НИ является классической задачей. Систематическое изложение методов их вычисления, существавших ранее, можно найти в [7, с. 314].

В [1] предложен весьма простой метод вычисления полиномов ${{a}_{n}}(x)$. Пусть

тогда где

Формулу (4) можно также использовать для вычисления НИ непосредственно по ее генератору $\omega (y)$ без громоздкого рекурсивного дифференцирования; или для быстрого вычисления генератора $\omega (y)$. Поскольку, согласно (3), имеем

то, подставив x = 0 в эту формулу, получим где величины ${{b}_{n}}$ вычисляются в (4); либо ${{b}_{n}}$ считаются известными, и тогда вычисляются величины ${{c}_{n}}$ одновременно с вычислением полиномов ${{a}_{n}}(x)$.

Приведем еще одну полезную формулу (см. [6, (84)]), связывающую НИ с ее генератором:

Как отмечалось выше, отображение (1) порождает ДДС ${{y}_{{n + 1}}} = F({{y}_{n}})$, $n \in {\text{Z}}$. При этом НИ ${{F}^{{(x)}}}(y)$, $x \in \mathbb{C}$, порождает НДС при условии, что формальные ряды ${{F}^{{(x)}}}(y)$ можно интерпретировать как некоторые функции. Если подразумевать под этим сходимость формальных рядов НИ, то это удается сделать крайне редко.

Насколько нам известно, все существующие примеры аналитических НИ даются формулой (см. [8])

соответствующей генератору $\omega (y) = {{y}^{{k + 1}}}$ (см. также [6, (72)]).

НИ дает пример динамической системы, инвариантные кривые которой можно сразу выписать в виде формальных рядов:

где $C = {{y}_{0}}$ – начальное значение (и первый интеграл), и $y \in \mathbb{R}$ или $y \in \mathbb{C}$. В последнем случае непрерывное “время” $x \in \mathbb{R}$ можно исключить из уравнений и рассматривать инвариантные кривые как функции $u = u({v})$ или ${v} = {v}(u)$, где $u = {\text{Re}}(y)$, ${v} = {\text{Im}}(y)$.

Например, рассмотрим простейшее из всех возможных нетривиальных отображений (1), F(y) = $y - {{y}^{2}}$, т.е. частный случай логистического отображения ($F(y) = y + A{{y}^{2}}$ приводится к этому виду перенормировкой):

Вычисления по формуле (4) дают генератор:

а также формальный ряд, представляющий инвариантные кривые НИ:

Напомним, что оба этих ряда расходятся (см. [1] и ссылки там).

Тем не менее можно легко проверить численно, что даже этот скромный отрезок формального ряда ${{(y - {{y}^{2}})}^{{( - x)}}}$ хорошо приближает инвариантную кривую $y = y(x)$ для $y(0) = C = 0.1$ по крайней мере в интервале $x \in [ - 2,1]$.

Для меньших значений $y(0) = C$, либо для большего отрезка ряда ${{(y - {{y}^{2}})}^{{( - x)}}}$ этот интервал можно значительно расширить в обе стороны (см. [1]).

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

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

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

В этой работе мы намереваемся сделать и то и другое.

3. ПЕРВЫЙ ИНТЕГРАЛ НЕПРЕРЫВНОЙ ИТЕРАЦИИ

К сожалению, применить теорему Бореля–Ритта напрямую к формальным рядам, соответствующим инвариантным кривым НИ, не представляется возможным, так как эти ряды зависят от двух переменных. Такие ряды ведут себя значительно сложнее обычных уже в случае их сходимости, а расходящиеся ряды от двух и более переменных вообще мало изучены.

Найдем функцию $H(x,y)$ (формальный ряд) такую, что ее линии уровня $H(x,y) = C$ будут соответствовать инвариантным кривым НИ. Положим формально

где “время” x соответствует начальному значению $y = y(x)$. Сдвиг по времени t будем считать малым.

Теорема 1. Для любого формального ряда (1) существует аналитический первый интеграл $H(x,y)$, линии уровня которого соответствуют инвариантным кривым непрерывной итерации при достаточно малых $y \in \mathbb{C}$.

Доказательство. Нужно доказать, что функция A(t) в (11) является плоской, т.е. имеет тождественный ноль в качестве степенного асимптотического разложения (в некотором секторе $\mathbb{C}$) для подходящей функции $H(x,y)$. Таким образом, первый интеграл (если существует) не является единственным. Однако он не может быть единственным даже в случае, если все ряды заведомо сходятся, поэтому неединственность $H(x,y)$ нас не должна беспокоить.

Запишем функцию $A(t)$ в (11) в виде

используя формулу (3), выражающую НИ с помощью ее генератора. Тогда

Поэтому в первом приближении необходимо

где G – это произвольная регулярная функция. Имея в виду неединственность первого интеграла, положим $G() = 1$.

Далее необходимо найти поправки к A(t) более высокого порядка по $t$. Однако в этом нет необходимости, так как, правильно выбрав пределы интегрирования в (12), мы получаем точное решение:

Доказательство этой формулы можно вывести из представления НИ (3) и формулы (6), однако это уже сделано ранее (см. [6, (88)]).

Интеграл (13) является пока формальным. Запишем его в форме, соответствующей обозначениям в (8), т.е. положим ${{y}_{0}} = C$ при времени x = 0 и $y = y(x)$ при времени x. Тогда получим другое выражение для инвариантной кривой (8), где переменные разделены:

Поскольку подынтегральное выражение в (14) является формальным рядом одной переменной (независимо от сходимости или расходимости ряда (1)), ему соответствует некоторая голоморфная функция в некотором секторе комплексной плоскости (теорема Бореля–Ритта, см. [9, с. 43]). Почленно интегрируя асимптотическое разложение $1{\text{/}}\omega (t)$ в (14) по пути $t \in [C,y]$, который следует выбрать в этом секторе, получим аналитическую функцию $x = x(y)$, а также утверждение теоремы. Теорема доказана.

Например, взяв генератор $\omega (y)$ (10) для логистического отображения (9), получим численно

причем тем точнее, чем меньше ${\text{|}}y{\text{|}}$, $y \in \mathbb{C}$. Дифференцируя эту формулу, получаем функциональное уравнение для генератора (10): откуда можно получить рекуррентную формулу для его коэффициентов.

Как мы отмечали, все существующие примеры аналитических НИ являются алгебраическими и даются формулой (7). Дадим пример трансцендентной НИ, соответствующей генератору $\omega (y) = {{y}^{2}} + {{y}^{3}}$. Формула (14) дает:

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

где ${{W}_{s}}()$ – это функция Ламберта, а $s$ – это ее ветвь: $s = 0$ при $y < 0$ и $s = - 1$ при $y > 0$.

Разложить формулу (16) в ряд по $y$ непосредственно не удается. Следует сперва разложить ее в ряд по $x$, а затем воспользоваться тождеством $y = W(yexp(y))$. Тогда получим ряд, соответствующий генератору по формуле (4).

Существует другой подход к вычислению первого интеграла $H(x,y)$, соответствующего инвариантным кривым отображения (1) (см. [1]). При этом мы не используем концепцию НИ, а рассматриваем отображение (1) как ДДС:

где переменная ${{x}_{n}}$ играет роль дискретного времени.

Очевидно, задача вычисления первого интеграла $H(x,y)$ такого, что $H({{x}_{{n + 1}}},{{y}_{{n + 1}}}) = H({{x}_{n}},{{y}_{n}})$, $n \in \mathbb{Z}$, эквивалентна вычислению НИ отображения (1), т.е. восстановлению НДС по ее дискретной редукции.

Сделаем преобразование масштабирования в системе (17):

где h – это малый параметр. Затем умножим первое уравнение (17) на h, а второе уравнение (17) разделим на h. Тогда получим ДДС:

Вычислим первый интеграл $H(x,y) = {\text{const}}$ для системы (18), а затем проделаем обратное масштабирование, что эквивалентно подстановке $h = 1$ в полученные ниже формулы.

Положим формально ${{x}_{n}} = x$, ${{y}_{n}} = y$ и рассмотрим тождество

Выбирая простейшее частное решение, получаем

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

Далее для простоты будем считать ${{c}_{1}} \ne 0$. В случае ${{c}_{1}} = 0$ масштабирование следовало делать по-другому. Тогда получим

Используем эту формулу в (19). Получим

Далее, как можно убедиться с помощью математической индукции, в разложении $H(x,y)$ по параметру h получим формальный ряд по y. Поэтому далее положим h = 1 и будем искать функцию $H(x,y)$ в виде

где

Тождество $H(x + 1,F(y)) = H(x,y)$, где $F(y)$ – это формальный ряд (1), дает функциональное уравнение для функции $U(y)$:

Предполагая, что функцию $U(y)$ можно каким-либо образом вычислить по уравнению (21), мы получили первый интеграл (20), т.е. восстановили НДС по ее дискретной редукции.

Однако на самом деле мы не получили ничего нового, кроме уравнения (21).

Предложение 1. Первый интеграл (20) совпадает с интегралом (14).

Доказательство. Положим $x = x(y)$ в формулах (14) и (20), затем продифференцируем их $y$ и приравняем результаты. Получим

Но это не что иное, как асимптотическое разложение подынтегрального выражения в (14), так как по формуле (4) получим

Поэтому если определить функцию U(y) как формальный ряд, полученный из ОДУ (22), т.е.

то получим требуемое утверждение. Предложение доказано.

Для общих отображений (1) масштабирование (т.е. искусственное введение малого параметра) следует подбирать индивидуально. Например, для итераций синуса ${{y}_{{n + 1}}} = sin{{y}_{n}}$, $n \in \mathbb{Z}$, следует взять

а затем умножить первое уравнение (17) на $h$, а второе уравнение (17) разделить на $\sqrt h $.

Однако проще сразу использовать асимптотическое разложение функции $1{\text{/}}\omega (y) = \ldots + U{\kern 1pt} '(y)$ в формуле (14), где $U(y)$ – это формальный ряд по натуральным степеням $y$.

Например, для итераций синуса оба подхода дадут результат:

где $U(y) = 79{\text{/}}3150{{y}^{2}} + 29{\text{/}}7875{{y}^{4}} + \ldots $ удовлетворяет функциональному уравнению

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

Функциональное уравнение (21) для функции U(y), введенное таким образом, дает, как мы увидим, эффективный способ ее вычисления (см. разд. 5).

4. АНАЛИТИЧЕСКИЕ СВОЙСТВА ПЕРВОГО ИНТЕГРАЛА НЕПРЕРЫВНОЙ ИТЕРАЦИИ

Как отмечалось в [10], существование аналитической функции, соответствующей асимптотическому разложению, представленному формальным степенным рядом, и вычисление этой функции – это две совершенно разные задачи.

Конструктивное вычисление аналитической функции $H(x,y)$, соответствующей НИ отображения (1), подразумевает также изучение аналитических свойств этой функции. В частности, изучение ее особенностей.

Ясно, что в общем случае это сделать невозможно, поэтому далее мы сосредоточимся на изучении аналитического интеграла $H(x,y)$ для функции $F(y) = y - {{y}^{2}}$, т.е. для логистического отображения (9).

Для этого отображения ${{c}_{1}} = - 1$, ${{c}_{n}} = 0$, $n > 1$. Поэтому первый интеграл (20) принимает вид

а функция U(y) удовлетворяет функциональному уравнению

Это те же уравнения, которые были получены в [1].

Напомним, что для коэффициентов формального ряда U(y) существует рекуррентное соотношение:

где где [x] – это целая часть x, и $C(n,m)$ – это биномиальный коэффициент.

Из результатов разд. 3 следует, что расходимость формального разложения (26) функции U(y) равносильна расходимости ряда (10) для генератора $\omega (y)$ НИ отображения (9). А расходимость ряда (10) влечет, в свою очередь, расходимость рядов, представляющих инвариантные кривые НИ этого отображения.

Доказательство расходимости рядов НИ для логистического отображения (9) являлось предметом диссертации (см. [11]). Однако мы не будем опираться на эти результаты, а докажем расходимость ряда для функции U(y) прямо из уравнения (25).

Теорема 2. Формальный ряд (26) всюду расходится, так как точка $y = 0$ является предельной точкой особенностей функции U(y).

Доказательство. Напомним, что функция U(y) является аналитической в некотором секторе, примыкающем к точке $y = 0$ (см. теорему 1).

Точка ${{y}_{0}} = 1$, очевидно, является особой для функции U(y), что следует из уравнения (25). По той же причине точки

являются особыми как прообразы точки y₀, т.е. $y - {{y}^{2}} = {{y}_{0}}$, $y = y_{1}^{ \pm }$.

Далее будем учитывать только точки в верхней полуплоскости ${\text{C}}$. Поэтому получим, что все точки

являются особыми.

Например, ${{y}_{{10}}} \approx - 0.12356942653420 + 0.037052550159283i$.

Но последовательность точек (27) стремится к нулю, так как имеет асимптотику

которая вычисляется методом неопределенных коэффициентов.

На самом деле, это в точности та же асимптотика, которую мы нашли в [1, (21)] с помощью первого интеграла (24). Нужно лишь заменить там t на –1/n и правильно подобрать константу C. Для последовательности (27), $C \approx - 0.46947082 + 2.50176142i$. Терема доказана.

Продолжим изучение особенностей функции U(y).

Рассуждая так же, как в доказательстве теоремы 2, получаем, что все точки в комплексной плоскости, которые за конечное число итераций $y \to y - {{y}^{2}}$ попадают в ноль, являются особыми для функции U(y).

Но это еще далеко не все ее особенности. Например, подставим в уравнение (25) $y = 1 + i$, затем $y = 1 - i$ и сложим эти уравнения. Тогда получим $0 = - 2 + i\pi $. Поэтому точки $y = 1 \pm i$ также являются особыми.

Рассуждая аналогично, получим, что все периодические орбиты отображения (9), $y = {{F}^{{(n)}}}(y)$, $n \in \mathbb{N}$, состоят из особых точек функции U(y). А следовательно, и все точки, которые за конечное число итераций $y \to y - {{y}^{2}}$ попадают в одну из точек орбиты, также являются особыми.

На самом деле множество всех особых точек функции U(y) – это не что иное, как множество Жюлиа $\mathcal{J}$ отображения (9) (см. [12, с. 53]), так как замыкание объединения указанных выше особых множеств также состоит из особенностей функции U(y).

На фиг. 1 показаны корни полиномов $y = {{(y - {{y}^{2}})}^{{(n)}}}$, для $n = 6,7,8$, и множество Жюлиа $\mathcal{J}$ для примерно 45 000 особых точек.

Если расположить корни ${{y}_{j}}$ полинома $y = {{(y - {{y}^{2}})}^{{(n)}}}$ в порядке возрастания аргумента ${\text{arg}}({{y}_{j}}$ – 1/2), то можно посчитать периметр ${{P}_{n}}$ многоугольника с данными вершинами. Оказалось, что ${{P}_{{n + 1}}}{\text{/}}{{P}_{n}} \approx 1.022$, $n > 6$. Поэтому периметр фигур на фиг. 1 возрастает неограниченно. А это означает, что множество Жюлиа $\mathcal{J}$ устроено фрактально, подобно звезде фон Коха. В частности, оно не является спрямляемой жордановой кривой.

Аналогичное множество Жюлиа для отображения

показано более подробно в [13, р. 42, фиг. 5а]. Там оно было названо “простой замкнутой кривой” и “довольно дикой жордановой кривой”. На наш взгляд, это множество также является фракталом.

Разумеется, точное вычисление фрактальной размерности множества $\mathcal{J}$ является отдельной и весьма сложной задачей.

Множество $\mathcal{J}$ является по определению замкнутым. Открытое множество $\mathbb{C}{{\backslash }}\mathcal{J}$ распадается на две отдельные компоненты: ${{\mathcal{F}}_{0}}$ и ${{\mathcal{F}}_{\infty }}$. Это множества Фату (см. [13, р. 39]), в которых любая точка $y$ попадает, соответственно, в ноль или бесконечность за бесконечное число итераций $y \to y - {{y}^{2}}$.

Множества Фату ${{\mathcal{F}}_{0}}$ и ${{\mathcal{F}}_{\infty }}$ – это, соответственно, внутренняя и внешняя части множества Жюлиа (см. фиг. 1).

По теореме Бореля–Ритта функция U(y) является голоморфной в некотором секторе, примыкающем к началу координат. Однако из уравнения (25) можно извлечь гораздо больше.

Теорема 3. Функция $U(y)$ является голоморфной в области Фату ${{\mathcal{F}}_{0}}$.

Доказательство. Рассмотрим любую точку ${{y}_{0}} \in {{\mathcal{F}}_{0}}$. По определению,

где ${{y}_{n}} = {{y}_{{n - 1}}} - y_{{n - 1}}^{2}$. Перепишем уравнение (25) в виде где логарифм мы упростили, так как он здесь однолистен. Затем просуммируем эти уравнения по $n \in {{\mathbb{N}}_{0}}$. Таким образом, имеем (пока формально)

В последней сумме существенно лишь то, что общий член имеет порядок $O(y_{n}^{2})$. Согласно [1, (21)], асимптотика последовательности {y_n}, $n \in \mathbb{N}$, дается формулой

где константа $C = C({{y}_{0}})$. Таким образом, общий член в сумме (29) имеет порядок $O(1{\text{/}}{{n}^{2}})$, поэтому этот ряд сходится для любой точки ${{y}_{0}} \in {{\mathcal{F}}_{0}}$.

Обозначим $y = {{y}_{0}}$ и перепишем формулу (29) в виде

где ${{G}_{n}}(y)$ – это голоморфная функция в некоторой окрестности точки $y$. Согласно данной выше оценке, этот ряд сходится равномерно. Согласно теореме Вейерштрасса (см. [14, р. 249]), функция U(y) также является голоморфной. Теорема доказана.

Формула (22) для отображения (9) имеет вид

где $\omega (y)$ дана в (10). Поэтому вычислим также функцию $U{\kern 1pt} '(y)$.

Ряд (31) (он же ряд (29)) можно почленно дифференцировать. Имеем

По правилу дифференцирования сложной функции получаем

Оценим асимптотику произведения в формуле (33). Согласно асимптотике (30), это произведение расходится к нулю. Точнее,

где мы оставили только главный член асимптотики, ${{y}_{k}} \asymp 1{\text{/}}k$. Поэтому ряд (33) также сходится (хотя мы знали это заранее по теореме Вейерштрасса, (см. [14, р. 249])).

Заметим, что мы только что еще раз доказали теорему 3, так как показали существование комплексной производной функции $U(y)$ в области Фату ${{\mathcal{F}}_{0}}$.

Формула (33) позволяет вычислить функцию $U{\kern 1pt} '(y)$ в некоторых точках в явном виде, так как ряд (33) обрывается при $y = 1{\text{/}}2$ и во всех ее прообразах. Например, $U{\kern 1pt} '(1{\text{/}}2) = 2$. Поскольку уравнение $y - {{y}^{2}} = 1{\text{/}}2$ имеет корни ${{y}_{ \pm }} = 1{\text{/}}2 \pm i{\text{/}}2$, то

Теорема 4. Генератор (10) $\omega (y)$ имеет счетное множество простых полюсов в области Фату ${{\mathcal{F}}_{0}}$ и не имеет там других особенностей. Эти полюса накапливаются вблизи множества Жюлиа $\mathcal{J}$. Кроме того, $\omega ({{y}_{ * }}) = 0$ на счетном всюду плотном множестве точек ${{y}_{ * }} \in \mathcal{J}$.

Доказательство. По формуле (32) получаем, что $\omega (1{\text{/}}2) = \infty $, т.е. точка $y = 1{\text{/}}2$ является особой для генератора (10). Поэтому все прообразы точки $y = 1{\text{/}}2$ также являются особыми для $\omega (y)$, что видно из функционального уравнения (15), которое можно записать в виде

Все прообразы точки $y = 1{\text{/}}2$ являются простыми полюсами по этой формуле. Кроме того, счетное множество этих полюсов функции $\omega (y)$ в области Фату ${{\mathcal{F}}_{0}}$ не может иметь точку сгущения внутри этой области. Иначе это была бы существенно особая точка для $\omega (y)$, а стало быть и для функции $U{\kern 1pt} '(y)$, которая, согласно теореме 3, является голоморфной в ${{\mathcal{F}}_{0}}$.

По той же причине $\omega (y)$ не может иметь других существенных особенностей либо точек ветвления внутри ${{\mathcal{F}}_{0}}$.

Других полюсов, кроме указанных, т.е. прообразов точки $y = 1{\text{/}}2$, у функции $\omega (y)$ внутри ${{\mathcal{F}}_{0}}$ быть не может. Иначе под действием итераций $y \to y - {{y}^{2}}$ эти полюса неизбежно попадут в любой сектор комплексной плоскости (с любым раствором), примыкающий к точке $y = 0$, что противоречит теореме Бореля–Ритта.

Это последнее утверждение следует из предложения 2, доказанного ниже.

Осталось найти нули функции $\omega (y)$.

Так как $\omega (0) = 0$ по определению, то из формулы (34) получаем, что все точки ${{y}_{ * }} \in \mathcal{J}$, которые за конечное число итераций попадают в ноль, также являются нулями $\omega (y)$.

Для периодических орбит в $\mathcal{J}$ система уравнений (34) имеет единственное решение $\omega (y) = 0$. Это же справедливо для прообразов этих точек.

Но эти множества точек всюду плотны в множестве Жюлиа $\mathcal{J}$. Терема доказана.

Разумеется, множество $\mathcal{J}$ также является особым для функции $\omega (y)$. Причем, по-видимому, функция $\omega (y)$ может принимать на нем любое комплексное значение.

Это видно из формулы (34): для каждой точки ${{y}_{ * }} \in \mathcal{J}$ выберем близкую к ней точку ${{y}_{0}} \in {{\mathcal{F}}_{0}}$. Имеем ${{y}_{n}} \to 0$ внутри ${{\mathcal{F}}_{0}}$, причем ${{y}_{n}} = O(1{\text{/}}n)$. Поэтому $\omega ({{y}_{n}}) = O(1{\text{/}}{{n}^{2}})$. Но произведение в (34) расходится к нулю с той же асимптотикой.

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

Тем не менее в следующем разделе мы покажем, что генератор $\omega (y)$, а также функции $U(y)$ и $U{\kern 1pt} '(y)$ эффективно вычисляются в своих регулярных точках внутри ${{\mathcal{F}}_{0}}$.

В завершение этого раздела приведем асимптотическую формулу для итераций отображения $y \to y - {{y}^{2}}$ в координатах $u = {\text{Re}}(y)$, $v = {\text{Im}}(y)$.

Для этого в формулу [1, (21)], т.е. в (30), нужно подставить $n = 1{\text{/}}t$, $C = A + iB$. Затем отделить вещественную и мнимую части разложения, что даст $u = t + O({{t}^{2}})$, ${v} = B{{t}^{2}} + O({{t}^{3}})$. Затем нужно обратить ряд для $u = u(t)$, что можно сделать методом неопределенных коэффициентов, либо применить алгоритм обращения степенно-логарифмических рядов, основанный на НИ (см. [1]). Получим степенно-логарифмический ряд $t = u + O({{u}^{2}})$. Затем этот последний нужно подставить в ряд для ${v} = {v}(t)$ и переразложить его по u. В результате получим

Непосредственно из интеграла (24) извлечь эту асимптотику не представляется возможным.

Другой способ вычислить эту асимптотику – это использовать первый интеграл (14). Запишем его в виде

и используем разложение (10). Затем подставим в него $y(x) = u(x) + i{v}(x)$ и отделим вещественную и мнимую части. Затем запишем (весьма громоздкое) уравнение $d{v}{\text{/}}du = \ldots $. Далее убеждаемся с помощью CAS, что это ОДУ имеет решение (35).

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

Предложение 2. Разложение (35) продолжается неограниченно и является степенным с одним вещественным параметром.

Доказательство. Мы дадим то, что называют “computer assisted proof”, так как без CAS воспроизвести это доказательство затруднительно.

Имеем, ${{y}_{n}} = {{u}_{n}} + i{{{v}}_{n}}$, ${{y}_{{n + 1}}} = {{u}_{{n + 1}}} + i{{{v}}_{{n + 1}}} = {{y}_{n}} - y_{n}^{2}$. Поэтому

(а): ${{u}_{{n + 1}}} = {{u}_{n}} - u_{n}^{2} + {v}_{n}^{2}$,

(б): ${{{v}}_{{n + 1}}} = {{{v}}_{n}} - 2{{u}_{n}}{{{v}}_{n}}$.

Мы уже имеем ${v} = B{{u}^{2}} + O({{u}^{3}})$, где B – это произвольная вещественная константа. Положим ${v} = B{{u}^{2}} + X{{u}^{3}} + O({{u}^{4}})$, где X – это неизвестный коэффициент. Тогда имеем

(в): ${{{v}}_{n}} = Bu_{n}^{2} + Xu_{n}^{3} + \ldots $,

(г): ${{{v}}_{{n + 1}}} - (Bu_{{n + 1}}^{2} + Xu_{{n + 1}}^{3} + \ldots ) = 0$.

Подставим (а) и (б) в уравнение (г), затем подставим туда (в) и ${{u}_{n}} = u$. Получим выражение, которое должно обращаться в ноль до любого порядка по u. Поэтому разложим его в ряд и найдем X из линейного уравнения при первой ненулевой степени $u$ (т.е. при u⁴). Получим X = B.

Далее положим ${v} = B{{u}^{2}} + B{{u}^{3}} + X{{u}^{4}} + O({{u}^{5}})$, и т.д.

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

Формула (35) описывает как стремление $y \to 0$ внутри области Фату ${{\mathcal{F}}_{0}}$ под действием итераций $y \to y - {{y}^{2}}$, так и асимптотику множества Жюлиа $\mathcal{J}$ вблизи нуля (для, возможно, канторова множества констант B).

5. КОНСТРУКТИВНОЕ ВЫЧИСЛЕНИЕ ПЕРВОГО ИНТЕГРАЛА НЕПРЕРЫВНОЙ ИТЕРАЦИИ

В предыдущем разделе мы показали, что функция U(y), а также и $U{\kern 1pt} '(y)$, и $\omega (y)$, однозначно вычисляются в своих регулярных точках в области Фату ${{\mathcal{F}}_{0}}$ с помощью сходящихся рядов. Область Фату ${{\mathcal{F}}_{\infty }}$ мы далее рассматривать не будем, так как там указанные функции будут заведомо многолистны.

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

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

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

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

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

Проще всего указанный алгоритм реализуется для функции $\omega (y)$ с помощью формулы (34). Для любой неособой точки ${{y}_{0}} \in {{\mathcal{F}}_{0}}$ вычисляются итерации (9), т.е. последовательность точек y_n до номера $n$, где величина |y_n| достаточно мала. Для точек, не слишком близких к границе области, количество итераций $y \to y - {{y}^{2}}$ невелико. Разумеется, следует предусмотреть остановку программы, если количество итераций велико (скажем, больше 1000). По мере вычисления точек ${{y}_{n}}$ накапливается произведение в (34). Затем вычисляется величина $\omega ({{y}_{n}})$ по разложению (10) и применяется (34).

Этим способом функция $\omega (y)$ (а значит, и $U{\kern 1pt} '(y)$, см. (32)) может быть вычислена (в принципе) с произвольной и контролируемой точностью.

На фиг. 2 показаны график функции $\omega (y)$ на вещественном интервале $y \in [0,1]$, и график функции ${\text{|}}\omega (1{\text{/}}2 + Rexp(i\phi )){\text{|}}$ на окружности вокруг особой точки $y = 1{\text{/}}2$ в комплексной плоскости при $R = 0.2$ и $R = 0.4$.

Тот факт, что функция $\omega (y)$ оказалась нечетной относительно своей особенности $y = 1{\text{/}}2$, виден из функционального уравнения (15), которое можно преобразовать к виду

Эта замена (${{y}_{n}} = {{z}_{n}} + 1{\text{/}}2$) приводит отображение (9) к виду ${{z}_{{n + 1}}} = - 1{\text{/}}4 - z_{n}^{2}$, что объясняет симметрию множества Жюлиа на фиг. 1.

Вычисления дают приближение функции $\omega (y)$ вблизи ее полюса $y = 1{\text{/}}2$:

Предложение 3. Функция $y(x)$ непрерывной итерации отображения (9) имеет в $\mathbb{C}$ счетное множество точек ветвления второго порядка.

Доказательство. Функция $\omega (y)$ имеет счетное множество простых полюсов (см. теорему 4). Запишем ОДУ (36) в локальных координатах в окрестности одного из ее полюсов

где $a$, $b$, $c$ – это некоторые комплексные константы. Тогда локально

Уравнение (36) может быть использовано для эффективного вычисления НИ логистического отображения для комплексных $y \in {{\mathcal{F}}_{0}}$ при некоторых ограничениях, связанных с предложением 3. Дело в том, что нас интересует главным образом вещественное “время” $x$. Поэтому кривая НИ должна избегать точек ветвления $y(x)$. Это возможно, например, при ${\text{Re}}({{y}_{0}}) < 1{\text{/}}2$, но невозможно при ${\text{Re}}({{y}_{0}}) > 1{\text{/}}2$ в связи с наличием полюса $y = 1{\text{/}}2$ у функции $\omega (y)$.

На фиг. 3 показан график НИ отображения (9) для комплексного начального значения ${{y}_{0}} = - 0.065 + i0.0135$, которое расположено вблизи границы области Фату ${{\mathcal{F}}_{0}}$, для вещественного времени $x \in [0,20]$. Показаны также дискретные итерации $y \to y - {{y}^{2}}$. Правый фиг. 3 – это фрагмент левого фиг. 3, на котором видно, что кривая НИ ведет себя контринтуитивно, демонстрируя осцилляции. Это, разумеется, вызвано близостью множества Жюлиа, где накапливаются особенности функции $\omega (y)$.

Если продолжить кривую НИ на фиг. 3 для бóльших значений x так, что ${\text{|}}y(x){\text{|}}$ станет достаточно мало, затем ввести новое локальное время x, т.е. начать с $x = 0$, то кривые НИ снова можно будет вычислять по их разложению, как мы это делали в [1] (см. также разд. 2).

Разумеется, аналитическая кривая на фиг. 3 вычисляется также из первого интеграла (24), но для этого нужно сначала вычислить функцию $U(y)$.

Следуя [1], введем новую функцию $V(y)$, которая не имеет особенностей на вещественной оси. Сделаем замену

в уравнении (25). Получим уравнение для функции $V(y)$:

Тогда интеграл (24) принимает вид

Следуя доказательству теоремы 3, запишем уравнение (38) в виде

и просуммируем его по $n$. Тогда получим

Функция $V(y)$ вычисляется по этой формуле внутри области Фату ${{\mathcal{F}}_{0}}$ точно так же, как описано в этом разделе для функции $\omega (y)$. Отличие состоит в том, что накапливается сумма, а не произведение, и $V(y)$ голоморфна внутри ${{\mathcal{F}}_{0}}$.

Функцию $V(y)$ можно также вычислить при $y < 0$. Запишем уравнение (38) в виде

и просуммируем его по $n$. Получим уравнение, аналогичное (40), по которому $V(y)$ вычисляется при $y < 0$ (см. (28) и (30)).

Функция $V(y)$ удобна тем, что она ограничена на интервале $y \in [0,1]$ и является четной по отношению к точке $y = 1{\text{/}}2$. Это видно после замены $y \to 1 - y$ в уравнении (38) (см. также (37)). Поэтому (38) можно записать в виде

Это уравнение позволяет вычислить асимптотику функции $V(y)$ при $y \to + \infty $ (см. [1]):

где $A \approx - 0.7641492$ – это однозначно определяемая константа.

Таким образом, функция $V(y)$ однозначно определена на всей вещественной оси (хотя в области Фату ${{\mathcal{F}}_{\infty }}$ эта функция будет многолистной).

На фиг. 4 показаны график функции $V(y)$, вычисленной по описанному алгоритму, а также график ее асимптотики (41), которую мы посчитали до члена $O(1{\text{/}}{{y}^{6}})$. Видно, что график асимптотики $V(y)$, $y \to + \infty $, практически совпадает с графиком $V(y)$, но не продолжается до малых $y$.

Асимптотика (41) была использована в [1] для вычисления асимптотики отображения (9) при ${\text{|}}y{\text{|}} \to \infty $.

В заключение этого раздела приведем значение функции $V(y)$ в точке ее глобального минимума при $y = 1{\text{/}}2$:

где, как мы полагаем, все знаки верны. Эти цифры совпадают для различных настроек алгоритма суммирования, а порядок погрешности (≈10^–32) соответствует учитываемой погрешности в асимптотическом разложении $V(y)$ в нуле.

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

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

Даже для простейшего из всех нетривиальных случаев отображения (1), т.е. для логистического отображения (9), формальный ряд генератора НИ, т.е. $\omega (y)$ (10) расходится. Мы это показали независимо от известных результатов (см. [11]) и весьма простым способом (см. теорему 2).

Тем не менее первый интеграл (14) определяет голоморфную функцию $x = x(y)$ во всей области Фату ${{\mathcal{F}}_{0}}$, при условии, что $C \in {{\mathcal{F}}_{0}}$. В частности, тождество

справедливо при любом ${{y}_{0}} \in {{\mathcal{F}}_{0}}$ и может быть проверено численно с любой (в принципе) точностью.

Таким образом, ДДС (9) является интегрируемой, хотя аналитический интеграл (39) выражается с помощью неизвестной ранее функции.

Кроме того, ДДС (9) дает пример интегрируемой системы, в которой реализуется динамический хаос. Это следует из хорошо известного хаотического поведения итераций дискретных отображений на их множествах Жюлиа (см. [13, c. 41]). В частности, любая точка ${{y}_{ * }} \in \mathcal{J}$ с вероятностью единица близка к некоторой точке периодической орбиты на $\mathcal{J}$ с неограниченным периодом.

Таким образом, интегрируемость в ДДС (а значит, и в НДС) не противоречит существованию хаоса в системе.

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

Дело в том, что нормализация системы всегда начинается с ее комплексификации. Но в данном случае “нормальная форма” ДДС (9) – это та же система, но в комплексных координатах.

Однако вещественная динамика в (9) совершенно тривиальна, в то время как комплексная динамика, как мы видели, совершенно нетривиальна.

И, наконец, система (9) дает еще один пример того, что так называемые “локальные методы” исследования особенностей применимы лишь к сравнительно простым случаям. В то время как исследование достаточно нетривиальной локальной задачи неизбежно приводит к необходимости изучать задачу глобально. См. также [10], где даны несколько аналогичных примеров.