Этот топик продолжает серию моих статей на Хабре, посвященных исследованию аттрактора Лоренца.

Часть 1. Критический взгляд на аттрактор Лоренца
Часть 2. Динамическая система Лоренца и вычислительный эксперимент
Часть 3. О существовании периодических решений в системе Лоренца
Часть 4. Три цикла в аттракторе Лоренца

Итак, рассмотрим нелинейную систему дифференциальных уравнений, введенную Эдвардом Лоренцом в 1963 году:

$$

где

$$

классические значения параметров системы.

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

В 70-х гг. 20-ого века Гукенхеймер, Уильямс и Йорке на основе результатов численных экспериментов сформулировали гипотезу о структуре аттрактора Лоренца при классических значениях параметров системы, однако, соответствие этой гипотезы структуре притягивающего множества системы (1) не было строго доказано. В 2000 году Стивен Смейл составил список из 18 наиболее значительных математических проблем 21-ого века. Проблема структуры аттрактора Лоренца была включена в этот перечень под номером 14. Уорвик Такер в своей работе 2002 года доказал, строго говоря, гиперболичность аттрактора в системе (1), т.е. аттрактор состоит из траекторий, всюду плотных на нем, вдоль которых близкие траектории экспоненциально разбегаются (континуум седловых траекторий); это и создает их хаотическое поведение. Тогда (как отмечает Д.В. Аносов в послесловии к книге Ж. Палиса и В. ди Мелу «Геометрическая теория динамических систем», Мир, 1986, с. 285) в аттракторе системы (1) "может (а не редко и должно) существовать бесконечное число асимптотически устойчивых периодических траекторий", но, на мой взгляд, их область притяжения может быть достаточно малой (трудно улавливаемой в численном эксперименте). На основе вычислений в работе В.С. Афраймовича, В.В. Быкова и Л.П. Шильникова «О возникновении и структуре аттрактора Лоренца» (ДАН СССР, 1977. — Т. 234, №2, с. 336-339) описана структура аттрактора — в нем также существует всюду плотное множество седловых циклов.

Как обычно отслеживают циклы в системе Лоренца? Используют символическую динамику — разбивают область в фазовом пространстве, содержащую аттрактор, на конечное число подобластей. Обозначая каждый элемент разбиения буквой, траектории на аттракторе, проходящие через соответствующие области, кодируются последовательностями таких символов. Если в последовательности имеется регулярность — повторяемость групп символов, — то мы имеем дело с циклом. Однако возращаемость траектории в некоторую окрестность своей части не говорит о ее замкнутости (часть 1). Вообще, критику результатов подобных вычислительных экспериментов можно найти в научной литературе (см., например, главу 4 «Can we trust in numerical computations of chaotic solutions of dynamical systems?», написанную французским математиком Рене Лози и опубликованную в 2013 году в книге «Topology and Dynamics of Chaos»).

В 2004 году Дивакар Вишванат опубликовал работу «The fractal property of the Lorenz attractor» (Physica D: Nonlinear Phenomena, 2004. — Vol. 190, iss. 1-2, с. 115-128), в которой привел начальные условия и периоды для трех циклов в аттракторе Лоренца. На мой взгляд, эта работа стала революционной в хаотической динамике, поскольку автор впервые показал конкретные значения фазовых координат и периода с достаточно большой точностью. Алгоритм вычислений был основан на методе Линдштедта-Пуанкаре (ЛП), на который (в отличие от методов численного интегрирования) не влияет устойчивость цикла, к которому строятся приближения.

Особенно впечатляет третий цикл.

Может сложится впечатление, что здесь несколько циклов, но нет — данный цикл имеет такую многоспиральную структуру из-за доминирования большого числа (порядка 100) гармоник.

Анализ работ Вишваната показал, что автор приводит описание алгоритма без программной реализации (на MATLAB, как написано в его статьях). При этом не ясно, как для ЛП-метода символьно решается получаемая неоднородная линейная система дифференциальных уравнений с периодическими коэффициентами (например, для уравнения Ван дер Поля это сделать можно без особых проблем). На письма с просьбой прислать исходные тексты программ автор не отвечает. Можно понять — он привел данные, которые можно проверить (часть 4), решая задачу Коши высокоточными численными методами (часть 2), а как конкретно все сделано — тут уже алгоритм закрыт, оставляя за собой право единственного авторства. Поиск работ по цитированиям показал, что потом никто не воспроизвел алгоритм Вишваната для отыскания периодических решений в системе Лоренца.

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

Метод рядов Фурье

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

$$

где h — заданное количество гармоник.

В силу правой части системы (1) составим невязки

$$

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

1. Продифференцировать по времени соответствующий тригонометрический полином;
2. Где имеются произведения фазовых координат, перемножить их, преобразовав произведения тригонометрических функций в суммы;
3. Привести подобные слагаемые для каждой функции cos() и sin() с соответствующим аргументом;
4. Поскольку мы работаем с количеством гармоник, равным h, отсечь от полученной невязки гармоники более высокого порядка;
5. Приравнять полученную невязку к нулю, т.е. коэффициенты при ее гармониках.

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

$$

а уравнений — на единицу меньше.

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

$$

параллельной плоскости xOy (сечение Пуанкаре). Таким образом, третья координата в начальном условии для искомых циклов равна величине r-1, откуда

$$

тогда дополнительное уравнение системы имеет вид:

$$

Других дополнительных сведений о периодических решениях в системе Лоренца я не встречал. Замечу, что для трех циклов, найденных Вишванатом, в начальном условии третьей координаты было взято число 27.

Далее приведен пример системы уравнений при h = 2:

$$

Беглый анализ для трех и более гармоник показал, что найти формулы общего вида некоторых уравнений системы пока затруднительно. Отметим, что для любого h такая система имеет решения

$$

соответствующие указанным положениям равновесия.

Пакет символьных вычислений Maxima

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

При разработке проекта мной был выбран математический пакет Maxima. Файл input_maxima.wxm с командами получения амплитуд и постоянных членов невязок при h = 2 представлен далее.

/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [wxMaxima: input start ] */
display2d:false$
x1:x10+c1c1*cos(1*omega*t)+s1c1*sin(1*omega*t)+
c1c2*cos(2*omega*t)+s1c2*sin(2*omega*t)$
x2:x20+c2c1*cos(1*omega*t)+s2c1*sin(1*omega*t)+
c2c2*cos(2*omega*t)+s2c2*sin(2*omega*t)$
x3:x30+c3c1*cos(1*omega*t)+s3c1*sin(1*omega*t)+
c3c2*cos(2*omega*t)+s3c2*sin(2*omega*t)$
assume(omega > 0)$
delta1:trigreduce(diff(x1,t)-(10*(x2-x1)),t)$
delta2:trigreduce(diff(x2,t)-(28*x1-x2-x1*x3),t)$
delta3:trigreduce(diff(x3,t)-(x1*x2-8/3*x3),t)$
expand(diff(delta1,cos(1*omega*t)));
expand(diff(delta1,sin(1*omega*t)));
expand(diff(delta1,cos(2*omega*t)));
expand(diff(delta1,sin(2*omega*t)));
expand(integrate(delta1,t,0,2*%pi/omega)*omega/(2*%pi));
expand(diff(delta2,cos(1*omega*t)));
expand(diff(delta2,sin(1*omega*t)));
expand(diff(delta2,cos(2*omega*t)));
expand(diff(delta2,sin(2*omega*t)));
expand(integrate(delta2,t,0,2*%pi/omega)*omega/(2*%pi));
expand(diff(delta3,cos(1*omega*t)));
expand(diff(delta3,sin(1*omega*t)));
expand(diff(delta3,cos(2*omega*t)));
expand(diff(delta3,sin(2*omega*t)));
expand(integrate(delta3,t,0,2*%pi/omega)*omega/(2*%pi));
/* [wxMaxima: input end ] */

Рассмотрим его подробнее. Комментарии в начале и конце нужны для того, чтобы файл можно было открыть в графическом фронтенде wxMaxima последних версий. Выражение display2d:false$ выключает многострочное рисование дробей, степеней и т.п. Знак $ позволяет вычислить результат выражения, но не выводить на экран (вместо ;). Функция trigreduce(выражение,t) свертывает все произведения тригонометрических функций относительно переменной t в комбинации сумм. Дифференцирование невязок по гармоническим функциям необходимо для получения соответствующих амплитуд. Функция expand(выражение) раскрывает скобки (выполняет умножение, возведение в степень, приводит подобные слагаемые). Для нахождения постоянных членов невязок применяется интегрирование на периоде, т.е. постоянный член k-ой невязки равен

$$

Чтобы при символьном интегрировании пакет не задавал вопрос о знаке частоты, дается команда

assume(omega > 0)$

Аналогично файл input_maxima.wxm формируется для любого количества гармоник. После выполнения данного скрипта пакет выведет в консоли символьные выражения для левой части системы алгебраических уравнений, которую я решал в нем же методом Ньютона. Например, для пяти гармоник скриптовый файл newton.wxm имеет следующий вид:


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

(%o3) [[omega = 3.984915779315225,x10 = -1.0530814406427574E-19,
x20 = -1.0530848010639959E-19,x30 = 23.17484126777741,
c1c1 = -5.73442859974501,s1c1 = 8.700967706873968,
c1c2 = -6.138073055934428E-19,s1c2 = -7.0677320406008846E-19,
c1c3 = 3.162874335852723,s1c3 = 2.061936729442591,
c1c4 = -8.5085912287645265E-19,s1c4 = 4.9192865298350146E-19,
c1c5 = 0.52386925730497,s1c5 = -0.691987635106,
c2c1 = -2.267166248701582,s2c1 = 10.98608920812201,
c2c2 = -1.210031439700999E-18,s2c2 = -2.9446761267244605E-19,
c2c3 = 5.627867598584345,s2c3 = -1.719199625236613,
c2c4 = -1.5456862362288386E-19,s2c4 = 1.8698116859813374E-18,
c2c5 = -0.8548869658075,s2c5 = -1.735775069972357,
c3c1 = -3.7683970563654851E-19,s3c1 = 5.9125298542126567E-19,
c3c2 = 7.278395231560698,s3c2 = -9.743072005078496,
c3c3 = 1.1948226901634459E-18,s3c3 = 1.5220232872406872E-18,
c3c4 = -3.45323649933811,s3c4 = -1.496599920543623,
c3c5 = 9.6317779065317723E-19,s3c5 = -1.0610022378452297E-18]]

Период, соответствующий вычисленной частоте, равен 1.576742309033. Полученное начальное условие

x(0) = -2.047685006587,
y(0) = 2.505814384075,
z(0) = 27.

Для 20-ти гармоник результаты следующие:


Здесь точность представления вещественного числа — 30 значащих цифр после точки, команда пакета

(fpprec:30,newtonepsilon:bfloat(10^(-fpprec+5)))$

Тогда полученный период T = 1.558652211165 (совпало 8 знаков после точки со значением, указанным Вишванатом), начальное условие

x(0) = -2.147375914135,
y(0) = 2.078143036771,
z(0) = 27.

Проверку этих чисел нужно делать численным решением задачи Коши (на больших отрезках времени перейти к высокоточным вычислениям). Я использую метод степенных рядов для построения неустойчивых решений системы Лоренца, описанный в части 2. Точность по ряду взял равной 1E-14, 80 бит под мантиссу вещественного числа. Приведу значения фазовых координат, вычисленных через период:

x(T) = -2.147466991199,
y(T) = 2.078421751536,
z(T) = 27.000633876853.

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

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

Автоматизация вычислений

Поскольку мы имеем дело с громоздкой системой уравнений, нужны программы, которые будут формировать скрипты для Maxima, запускать пакет и читать результаты вычислений. Такие программы были написаны на языке C++. Все исходные тексты, а также скрипты для пакета, собраны в проект periodic_sols на GitHub. Разберем их.

Взаимодействие с пакетом будем вести через перенаправление ввода/вывода на файлы. Пример команды запуска пакета:

maxima < input_maxima.wxm > output_maxima.txt

Для чтения из текстового файла результатов выполнения команд имеется функция ReadFromMaxima() (файлы maxima.h и maxima.cpp).

void ReadFromMaxima(std::fstream &f, std::string &res)
{
    std::string s;
    res = "";
    while(1)
    {
        f >> s;
        if(s[1] == '%' && s[2] == 'o')
        {
            while(1)
            {
                f >> s;
                if(s[1] == '%' && s[2] == 'i')
                    break;
                res += s;
            }
            break;
        }
    }
}

Первый аргумент — поток, связанный, например, с файлом output_maxima.txt, второй — строка, в которую считается результат. Если он многострочный (Maxima разбивает длинные строки на части), то результирующие строки склеятся в одну строку. Для запуска пакета применяется стандартная библиотечная функция system() из cstdlib.

Проект periodic_sols состоит из трех программ: periodic_sols.cpp, init_newton.cpp и calc_init_point.cpp. Для их сборки запускаем

./compile

Перед запуском ./periodic_sols необходимо подготовить файлы input.txt и system.txt. Формат первого файла (я не стал заморачиваться с интерфейсом — все программы консольные):

1-ая строка — порядок системы (в нашем случае он равен 3, в программе хранится в переменной n, это сделано не только для рассматриваемой системы);
2-ая строка — количество гармоник;
3-я строка — точность представления вещественного числа;
4-ая строка — номер фазовой координаты для дополнительного уравнения (в нашем случае он равен 3);
5-ая строка — значение этой координаты (в нашем случае оно равно 27);
6-ая строка — начальное приближение по частоте;
7-ая строка — начальное приближение по постоянным членам;
8-ая строка — начальное приближение по амплитудам косинусов (по синусам берется 0).

В файле system.txt находится правая часть динамической системы:

10*(x2-x1)
28*x1-x2-x1*x3
x1*x2-8/3*x3

В результате получается файл newton.wxm с командами для метода Ньютона. Этот файл подается на вход пакету. Программа init_newton.cpp создает файл new_newton.wxm на базе файла newton.wxm, записывая туда результаты вычислений для предыдущего количества гармоник (переменная prev_harmonics), меньшего или равного текущему числу (переменная harmonics). Это необходимо, когда мы хотим улучшить приближение к периодическому решению. Здесь алгоритм такой:

1. Решив систему нелинейных уравнений для малого количества гармоник (prev_harmonics), формируем файл res_newton.txt по результатам вычислений;
2. В файле input.txt увеличиваем количество гармоник и запускаем ./periodic_sols для формирования новой системы;
3. Запускаем программу ./init_newton, которая переносит вычисленные значения на соответствующие места начального приближения для новой системы, а недостающие заполняет нулями.

Программа calc_init_point.cpp на базе файла res_newton.txt вычисляет начальные условия и период, записывая в файл init_point.txt. Буква «b» означает то же самое, что и буква «E» в представлении вещественного числа.

Что дальше?

1. Самая затратная по времени операция — символьное интегрирование. Например, для 120 гармоник время формирования системы — более 2-х суток. Тут можно распараллелить вычислительный процесс на три нода, но это значительного эффекта не даст (Вишванат, вообще, вычислял 1024 гармоники). Здесь нужно изменить принцип формирования постоянных членов. Например, вычесть из общего выражения полученные до этого выражения для каждой гармоники. Кстати, при решении системы нелинейных уравнений методом Ньютона матрица Якоби не обращается, в пакете используется LU-разложение для решения системы линейных уравнений на каждой итерации метода;
2. Сейчас я ищу аналогии в получаемых системах нелинейных уравнений для разных количеств гармоник, чтобы система формировалась сразу. Если кто-нибудь тоже хочет попробовать, то здесь имеется архив с системами для 2-7 гармоник. Если нужно еще, пишите в комментариях, добавлю;
3. В директорию T__undefined помещены приближения к пока еще неизвестному периодическому решению системы Лоренца (номерные директории внутри соответствуют найденному количеству гармоник). Пока еще недосчитал. Если метод для данного случая все-таки сходится, то период цикла будет больше периода третьего цикла Вишваната;
4. Как писал выше, главная проблема метода Ньютона — это выбор начального приближения. Если будет получен общий вид системы, то хорошо бы отделить от нее решения, соответствующие положениям равновесия, т.е. получить новую систему, решения которой — только искомые периодические решения. Тогда начальное приближение гораздо проще подбирать. Отмечу, конце исходного файла periodic_sols.cpp приведен комментарий с первоначально подобранными начальными приближениями для первого цикла Вишваната (как потом выяснилось) и неизвестного периодического решения.