Динамическая система Лоренца и вычислительный эксперимент

    image

    Данный пост является продолжением моей статьи [1] на Хабрахабре об аттракторе Лоренца. Здесь рассмотрим метод построения приближенных решений соответствующей системы, уделив внимание программной реализации.

    Динамической системой Лоренца является автономная система обыкновенных дифференциальных уравнений третьего порядка

    image

    где image, r и b являются положительными числами. Для исследования поведения решений системы обычно берут классические значения параметров системы, т.е.

    image

    Как было отмечено в статье [1], в этом случае в системе (1) имеет место неустойчивость ее решений на аттракторе. По сути, это делает некорректным применение классических численных методов на больших отрезках времени (а на таких отрезках и строятся предельные множества динамических систем). Одним из вариантов преодоления этой проблемы является переход к высокоточным вычислениям, но такой подход ставит исследователя в жесткие рамки: во-первых, малая степень свободы для уменьшения ошибки (изменение величины шага image интегрирования и точности представления вещественного числа для управления вычислительным процессом), во-вторых, большой объем вычислений при очень малых image.

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

    image

    как

    image

    где коэффициенты image, image и image определяются начальными условиями; также получена длина (шаг) image отрезка сходимости рядов (2):

    image

    если image, то image, иначе image,

    image — любое положительное число.

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

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

    По формулам (3) вычисляем image, image и image до такого значения i, при котором имеет место оценка

    image

    Понятно, что значительного накопления ошибки интегрирования на длинных временных отрезках нам и здесь не избежать, поэтому будем реализовывать высокоточные вычисления с плавающей точкой на базе библиотеки GNU MPFR Library, а точнее, C++-интерфейса библиотеки MPFR для работы с вещественными числами произвольной точности — MPFR C++. Она удобна тем, что в ней имеется класс mpreal с перегруженными арифметическими операциями и дружественными математическими функциями. Для установки библиотеки в Ubuntu Linux выполним

    sudo apt-get install libmpfrc++-dev
    

    Кроме пакета libmpfr-dev менеджеров пакетов потянет еще и libgmp-dev. Это devel-пакет библиотеки GNU Multiple Precision Arithmetic Library или GMP (MPFR является ее расширением).

    Рассмотрим пример кода на C++ вычисления значений фазовых координат в конечный момент времени, а также проверки найденных значений.

    #include <iostream>
    #include <vector>
    using namespace std;
    
    #include <mpreal.h>
    using namespace mpfr;
    
    // Точность по степенному ряду
    #define eps_c   "1e-50"
    
    // Параметры системы уравнений Лоренца
    #define sigma   10
    #define r       28
    #define b       8/(mpreal)3
    
    // Количество бит под мантиссу вещественного числа
    #define prec    180
    
    // Как считать шаг: 1 - по оценке отрезка сходимости, 0 - заданная величина
    #define FL_CALC 1
    
    // Фиксированный шаг по времени
    #define step_t  "0.02"
    
    // Функция возвращает длину отрезка сходимости степенного ряда
    mpreal get_delta_t(mpreal &alpha0, mpreal &beta0, mpreal &gamma0)
    {
        mpreal h2 = (mpfr::max)((mpfr::max)(fabs(alpha0), fabs(beta0)), fabs(gamma0)),
               h1 = (mpfr::max)((mpfr::max)(2*sigma, r+2*h2+1), b+2*h2+1);
        mpreal h3 = h2 >= 1 ? h1 * h2 : (mpfr::max)((mpfr::max)(2*sigma, r+2), b+1);
        return 1/(h3 + "1e-10");
    }
    
    // Функция вычисления значений фазовых координат в конечный момент времени
    // x, y и z - координаты начальной точки; T - длина отрезка интегрирования;
    // way - направление поиска решений: 1 - вперед по времени, -1 - назад по времени
    void calc(mpreal &x, mpreal &y, mpreal &z, mpreal T, int way = 1)
    {
        cout << "\nВ начальный момент времени:\n\ndx/dt = " << sigma*(y-x) << "\ndy/dt = " <<
                r*x-y-x*z << "\ndz/dt = " << x*y-b*z << endl;
    
        mpreal t = 0, delta_t, L, p, s1, s2;
        bool fl_rp;
        do
        {
            if(FL_CALC)
                delta_t = get_delta_t(x, y, z);
            else
                delta_t = step_t;
            t += delta_t;
            if(t < T)
                fl_rp = true;
            else if(t > T)
            {
                delta_t -= t-T;
                fl_rp = false;
            }
            else
                fl_rp = false;
    
            vector<mpreal> alpha, beta, gamma;
            alpha.push_back(x);
            beta.push_back(y);
            gamma.push_back(z);
    
            int i = 0;
            L = sqrt(alpha[0]*alpha[0] + beta[0]*beta[0] + gamma[0]*gamma[0]);
            p = way * delta_t;
            while(L > eps_c)
            {
                // Вычисляем новые коэффициенты степенных рядов
                s1 = s2 = 0;
                for(int j = 0; j <= i; j++)
                {
                    s1 += alpha[j] * gamma[i-j];
                    s2 += alpha[j] * beta[i-j];
                }
                alpha.push_back(sigma*(beta[i] - alpha[i])/(i+1));
                beta.push_back((r*alpha[i] - beta[i] - s1)/(i+1));
                gamma.push_back((s2 - b*gamma[i])/(i+1));
    
                i++;
    
                x += alpha[i] * p;
                y += beta[i] * p;
                z += gamma[i] * p;
                L = fabs(p) * sqrt(alpha[i]*alpha[i] + beta[i]*beta[i] + gamma[i]*gamma[i]);
                p *= way * delta_t;
            }
        }
        while(fl_rp);
    
        cout << "\nКоординаты в конечный момент времени:\nx = " << x.toString() << "\ny = " <<
                y.toString() << "\nz = " << z.toString() << endl;
        cout << "\nЗначения производных:\n\ndx/dt = " << sigma*(y-x) << "\ndy/dt = " <<
                r*x-y-x*z << "\ndz/dt = " << x*y-b*z << endl;
    }
    
    int main()
    {
        mpreal::set_default_prec(prec);
        cout << "Машинный эпсилон = " << machine_epsilon() << endl;
    
        mpreal T;
        cout << "\nВведите длину отрезка времени > ";
        cin >> T;
    
        mpreal x, y, z;
        cout << "\nx0 > ";
        cin >> x;
    
        cout << "y0 > ";
        cin >> y;
    
        cout << "z0 > ";
        cin >> z;
        cout << endl;
    
        calc(x, y, z, T);
        cout << "\n\n*** Проход назад ***\n";
        calc(x, y, z, T, -1);
    
        return 0;
    }
    

    Для компиляции файла lorenz.cpp с этим кодом выполним следующее

    g++ lorenz.cpp -std=gnu++11 -lmpfr -lgmp -o lorenz
    

    Как видно из листинга программы, для хранения значений коэффициентов степенных рядов используются векторы из библиотеки STL. Было принято, что количество бит под мантиссу вещественного числа равно 180. Диапазон изменения экспоненты — по умолчанию от image до image (на 32-битных платформах). Функция machine_epsilon() определения машинного эпсилон image была использована для того, чтобы узнать его значение для заданного представления вещественного числа. В данном случае

    image

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

    image

    Отрезок времени, на котором производилось вычисления, — image. Первоначально точность оценки общего члена ряда бралась равной image. Этого достаточно, чтобы получить результат

    image;

    производные:

    image

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

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

    Указанная разрядность для вещественного числа была выбрана с целью отследить не только возврат траектории системы Лоренца в окрестность начальных условий [1, рис. 1], но и для прохода назад по времени от конечной точки к начальной по дуге траектории (равенство параметра way функции calc() -1). Тогда в расчетах нужно брать image. В результате получаются значения, совпадающие с начальными до седьмого знака после запятой (метод toString() класса mpreal, вызванный без параметров, говорит о преобразовании числа в строку с максимальным количеством знаков после запятой):

    x = 13.412656286837273085165416945301946328440634370684244
    y = 13.4643000297481126631507883918720904312092673686014399
    z = 33.4615641630148784946354299167181879731357599130041067

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

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

    В работе [5, с. 90, 91] для исследования траекторий системы (1) применяется метод Эйлера с переменным шагом image интегрирования. Величина image выбирается с отслеживанием ошибки на текущем шаге (т.е. локальный контроль), используя интервальную арифметику. Однако нет проверки общей ошибки интегрирования. Проход назад по времени, описанный выше, нам гарантирует правильность построения приближенного решения. Варьируя точностью оценки общего члена ряда, можно уменьшить ошибку на шаге, что не позволяет метод Эйлера. Также в рассмотренной модификации метода степенных рядов преимуществом перед общей схемой метода рядов Тейлора является быстрый расчет по формулам (3) коэффициентов разложения по сравнению с процедурой символьного дифференцирования правых частей уравнений системы (кроме того, в нелинейном случае под символьные выражения требуется много памяти для их хранения при вычислении производных высших порядков).

    Таким образом, зная состояние системы Лоренца в прошлом, мы с достаточной степенью точности можем предсказать поведение ее траекторий в течение длительных интервалов времени, а также вернуться назад. По сути здесь нарушается формальное определение хаоса [6, с. 118, 119].

    Литература


    1. Пчелинцев А.Н. Критический взгляд на аттрактор Лоренца, 2013. Хабрахабр.
    2. Pchelintsev A.N. Numerical and physical modeling of the dynamics of the Lorenz system // Numerical Analysis and Applications, 7(2): 159-167, 2014. DOI: 10.1134/S1995423914020098.
    3. Магницкий Н.А., Сидоров С.В. Новые методы хаотической динамики. — М.: Едиториал УРСС, 2004.
    4. Пчелинцев А.Н. О моделировании динамики системы Лоренца // Международная конференция «Колмогоровские чтения VI. Общие проблемы управления и их приложения (ОПУ-2013)», Тамбовский государственный университет им. Г.Р. Державина (Тамбов, 2013 год). Math-Net.Ru.
    5. Tucker W. A Rigorous ODE Solver and Smale's 14th Problem // Foundations of Computational Mathematics, 2(1): 53-117, 2002.
    6. Берже П., Помо И., Видаль К. Порядок в хаосе. О детерминистком подходе к турбулентности. — М.: Мир, 1991.
    Share post

    Similar posts

    Comments 20

      +6
      Некоторые придирки.

      Почему это формальное определение хаоса нарушается? Система детерминированная, значит и странный аттрактор — детерминированный хаос (вот например есть лекции, поновее чем [6], В.С. Анищенко, «Знакомство с нелинейной динамикой», 2002). Ничего необыкновенного в том, что мы можем просмотреть траектории даже в аттракторе я не вижу.

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

      ничего не значит, все дело как раз в «качественной» теории. То есть она предсказывает, а численный эксперимент может как-то подтверждать, да и то — условно.
        0
        Я сослался на работу, ставшей классической. В.С. Анищенко в этой книге называет детерминированным хаосом то, когда траектория не замыкается и воспроизводит некий апериодический процесс (стр. 64). Получается, что, например, такой классический объект теории дифференциальных уравнений, как почти периодическая функция, является хаотической? :) А в системе Лоренца предельные решения описываются рекуррентными траекториями, при этом почти периодическое решение (как и периодическое решение) — их частный случай.

        ничего не значит, все дело как раз в «качественной» теории. То есть она предсказывает, а численный эксперимент может как-то подтверждать

        Да, это и имелось в виду, сейчас поправлю, спасибо за замечание.
        0
        Интересно, в свое время игрался с AUTO и открыл для себя то, что странный аттактор Лоренца по сути дела является набором неустойчивых циклов, между которыми «пляшет» траектория системы. Тогда же и возникла мысль, что если проводить численное интегрирование с абсолютной точностью, то никакого странного аттрактора бы и не было.
          0
          Если стартовать из окрестности точки (0,0,0), то рано или поздно траектория всё равно в неё вернётся — чтобы сделать новую петлю в том же полупространстве, откуда пришло, или обойти стороной сепаратрису, утыкающуюся строго в начало координат, и уйти во второе полупространство. Ничего необычного в этом нет, равно как и нет разрушения существующего детерминированного хаоса.
          Тот факт, что системе Лоренца можно сопоставить одномерное отображение, подобное известному отображению «тент», если взять последовательность максимумов переменной z, занумеровать её и построить зависимость z[n+1] ( z[n] ). Получается картина, качественно и, если отмасштабировать, то и количественно подобная отображению:
          z[n+1] = z[n]/2 при z[n] < 1/2, (1 — z[n])/2 при z[n] > 1/2.
          Получается вполне определённая последовательность, в которой все элементы до бесконечности определяются начальным состоянием. А вот то, хаотичное будет поведение или нет, зависит, конечно, от самого начального состояния. Не помню, как с «тентом», но с более простым отображением было просто достаточно задать начальное положение иррациональным числом, и система никогда бы уже в него не вернулась. В то время как при рациональном значении получился бы цикл, определяемый только периодичностью двоичного представления этого рационального числа.
          Да, известно так же, что в системе Лоренца есть циклы, но они неустойчивы, а суметь попасть на них крайне трудно — множество этих циклов так же фрактально, как и всем известная «бабочка» хаотической траектории (Кузнецов в своих книгах называл это множество «странным репеллером», т.к. оно в силу неустойчивости «отталкивает» траекторию). Но попасть в них всё-таки возможно.
          А говорить о нарушении хаоса как минимум преждевременно.
            0
            Да, известно так же, что в системе Лоренца есть циклы

            Строго существование циклов в системе Лоренца при классических значениях ее параметров до сих пор не доказано. Здесь нужен специальный математический аппарат. На сегодняшний день удалось показать лишь существование предельного цикла при больших значениях параметра r методом усреднения. Как говорил один знакомый математик, «если считается, что в системе Лоренца есть хотя бы один цикл, то укажите мне его начальную точку и период». Все заключения, которые есть в литературе о периодических решениях в системе, основаны на гипотезе, возникшей из численного эксперимента.
            0
            Динамической системой Лоренца является автономная система обыкновенных дифференциальных уравнений третьего порядка

            Судя по представленным ниже уравнениям, ОДУ в системе первого порядка. Т.е. это система трех диф. уравнений первого порядка.
              +3
              Вы говорите о порядке системы относительно каждой из функций x(t), y(t) и z(t). Число 3 является порядком системы — см. стр. 26 классической книги Л.С. Понтрягина «Обыкновенные дифференциальные уравнения».
                +1
                Тогда фразу следует перестроить, т.к. «система обыкновенных дифференциальных уравнений третьего порядка» имеет неопределенность к чему относится указание порядка (к дифф. уравнениям или же к системе).
                Пример: «система третьего порядка обыкновенных дифференциальных уравнений» — пусть и корявая, но уже не требует разбора контекста.
                  +1
                  Дифурщики так не говорят.
                    0
                    Лично я стараюсь к миниму свести неоднозначности в высказываниях (не всегда получается).
                    Про корявость фразы я специально написал. Быстро устранить неоднозначность в данном случае не мог. Пришлось лишь указать на оную.
                    Касаемо «диффурщиков», могу сказать за себя. Системы ДУ занимают большУю часть времени в моих мыслях. Всегда стараюсь разделять порядок системы (или, например, порядок модели Пространства Состояний) и порядок самих уравнений в системе.
                    На моей практике чаще всего приходится рассматривать одно ДУ высокого порядка, которое раскладывается в Форму Коши, т.е. в систему ДУ первого порядка. В этом случае порядок системы станет равен порядку исходного ДУ высокого порядка. Ожнако, даже в этом случае стараюсь не допускать, чтобы эти два «порядка» в смысле фразы смешивались.
              0
              во-первых, малая степень свободы для уменьшения ошибки (изменение величины шага image интегрирования и точности представления вещественного числа для управления вычислительным процессом)

              Некорректное, на мой взгляд, использование термина "степень свободы".
              Даже в мат. статистике этот термин означает скорее количество независимых случайных переменных в некоем случайном векторе. Здесь же оно использовано как характеристика числовой обусловленности, если я правильно понял. Т.е. допустимый диапазон варьирования кванта времени численного интегрирования, в котором сохраняется устойчивость решений.
                0
                Здесь имеется в виду количество независимых переменных, которые управляют процессом численного интегрирования системы: разрядность вещественного числа, величина шага, максимальный порядок производной при усечении степенного ряда.
                  0
                  Понял. Спасибо за разъяснения.
                0
                поэтому будем реализовывать высокоточные вычисления с плавающей точкой на базе библиотеки GNU MPFR Library

                Реализовать «вычисления с плавающей точкой» невозможно на библиотеке вещественных типов с фиксированной точкой. Видимо имелась в виду возможность изменения точности математических операций, увеличивая ее, когда нужно использованием представлений чисел с бОльшим числом разрядов дробной части.
                  0
                  Читайте внимательно описание библиотеки MPFR www.mpfr.org
                    0
                    The MPFR library is a C library for multiple-precision floating-point computations with correct rounding.

                    ХМММ Видимо я по-другому понимаю термин «floating-point» в данном случае. Для меня он означает, что конкретное число в форме floating-point задается неким дробным числом помноженным на степенную функцию. Производя действия над ним, мы меняем показатель степенной функции и тем самым смещаем разделитель дробной и целой частей.
                    Здесь же, если я правильно понимаю, позиция разделителя задается в момент создания экземпляра переменной и не меняется на протяжении жизни. Если результат действия может выйти за пределы заданного при создании переменной диапазона, то под хранение результата создается переменная увеличенной разрядности либо уменьшенной точности.
                    С какой-то точки зрения сами операций, наверное, можно обозначить термином Floating-Point, но я бы обозначил его Variable-Precision Operations.
                      0
                      Разрядность не меняется, просто точка в записи числа может быть помещена где угодно относительно его цифр. Поэтому и плавающая точка.
                  +1
                  Я внес небольшие изменения «косметического» характера в код. Спасибо за замечания Павлу Голобородько — разработчику MPFR C++.
                    0
                    Код и текст топика я опять поправил — в MPFR C++ есть встроенная функция определения машинного эпсилон, она точнее его вычисляет; оставил для экспоненты диапазон по умолчанию (он шире). И метод toString() можно вызывать без параметров.
                      0
                      Тут были ко мне вопросы о несовпадении координат точки в конечный момент времени. Предыдущая версия моей программы производила вычисления не совсем до конечного момента времени, а до T+delta_t (немножко пересчитывала), поэтому результаты сходились до первого знака после запятой с результатами, полученными в других программах. Я немного поправил код, параметры и результаты вычислений. Принципиально ничего не изменилось.

                      Only users with full accounts can post comments. Log in, please.