CC BY
22Научно-образовательный журнал для студентов и преподавателей «StudNet» №1/2021
ИССЛЕДОВАНИЕ МОДЕЛЕЙ ДИНАМИЧЕСКОГО ХАОСА НА ПРИМЕРЕ МОДЕЛИ ЛОРЕНЦА
STUDY OF MODELS OF DYNAMIC CHAOS ON THE EXAMPLE OF THE
LORENTZ MODEL
УДК 004.5
Макаров Д.А., студент, [email protected], Россия, 105005, г. Москва, МГТУ им. Н.Э. Баумана, кафедра «Системы обработки информации и управления»
Шибанова А.Д., студент, эЛа [email protected], Россия, 105005, г. Москва, МГТУ им. Н.Э. Баумана, кафедра «Системы автоматизированного проектирования»
Makarov D.A., student, [email protected], Russia, 105005, Moscow, MSTU N.E. Bauman, Department of Information Processing and Management Systems
Shibanova A.D., student, [email protected], Russia, 105005, Moscow, MSTU N.E. Bauman, Department of Computer-Aided Design Systems
Аннотация
В данной статье была исследована модель Лоренца (одна из наиболее известных моделей динамического хаоса) и построена классификация ее режимов с помощью численных методов решения задачи Коши (в частности -метода Эйлера, неявного метода Эйлера, а также метода Адамса-Башфорта-Моултона 4-го порядка) и случайно сгенерированных начальных условий. Были сделаны выводы о временных затратах, возникающих при решении
задачи каждым из методов, а также о том, как меняется фазовый портрет в зависимости от шага по времени й.
Summary
In this article, the Lorentz model (one of the most famous models of dynamic chaos) was investigated and a classification of its regimes was constructed using numerical methods for solving the Cauchy problem (in particular, the Euler method, the implicit Euler method, and the Adams-Bashfort-Moulton method of the 4th order) and randomly generated initial conditions.
Conclusions were made about the time costs arising in solving the problem by each of the methods, as well as about how the phase portrait changes depending on the time step й.
Ключевые слова: динамическая система, задача Коши, метод Эйлера, численный метод, аттрактор Лоренца.
Key words: dynamical system, Cauchy problem, Euler's method, numerical method, Lorentz attractor.
ОПРЕДЕЛЕНИЕ СТАЦИОНАРНЫХ ПОЗИЦИЙ ЗАДАННОЙ
СИСТЕМЫ ОДУ
Дана система ОДУ 1-го порядка: ^ = f(x), х — [x(t),y(t),z(t)]T
d 'х'
dt У —
-Z.
а(у — х) x(r — z) — у ху — bz
(1)
где а = 10 и b = 8/3.
Стационарной позицией динамической системы является постоянная во времени траектория x*(t) = const. Очевидно, что для стационарной позиции dx/dt = 0, из чего следует, что стационарные позиции можно найти, решив нелинейное в общем случае уравнение f(x) = 0. Тогда можно получить следующую систему уравнений:
о(у — х) = 0
- х(г-г)-у = 0 (2)
ху — Ьг = 0
Заметим, что симметрия уравнений Лоренца, заключающая в том, что вид этих уравнений не изменится, если одновременно сменить знаки у х и у, означает, что любое образование в фазовом пространстве либо обладает той же симметрией (то есть превращается само в себя при замене х —» —х и у —» —у), либо имеет такое же образование в качестве симметричной копии. Из 1-го уравнения системы (2) получаем: у = х.
Из 2-го уравнения системы (2) получаем:
у(г — г) = у
(3)
Уравнение (4) позволяет получить частное решение системы (2): X = у = 2 = 0
(4)
Выразим х из уравнения (4): 2 = г — 1
(5)
Тогда 3-е уравнение системы (3) можно записать так: у 2 — Ь(г — 1) = 0
(6)
Определим дискриминант квадратного уравнения (6): Б = —4Ь(г — 1)
(7)
Тогда решением уравнения (7) будет:
У = =)/ь{г-1)
(8)
Таким образом, получаем три возможных стационарных точки для заданной системы ОДУ:
1. х = у = 2 = 0;
2. х = — 1), у = — 1), г = г — 1
3. х = —^Ь(г — 1), у = —^Ь(г — 1), г = г — 1
Можно сделать вывод, что все возможные траектории системы ОДУ будут приближаться к этим трем точкам.
РЕАЛИЗАЦИЯ ФУНКЦИЙ ЧИСЛЕННЫХ МЕТОДОВ РЕШЕНИЯ
ЗАДАЧИ КОШИ
Чтобы построить траектории заданной (1) системы ОДУ, воспользуемся численными методами решения задачи Коши: методом Эйлера, неявным методом Эйлера, а также методом Адамса-Башфорта-Моултона. Сначала рассмотрим самый простой (и, к сожалению, не самый точный) метод решения систем ОДУ - метод Эйлера. Это явный, одношаговый метод 1-го порядка точности. Он основан на аппроксимации интегральной кривой кусочно- линейной функцией (ломаной Эйлера). Запишем формулировку метода Эйлера [1]: = а,
= + к/(^, I = 0,1 ... т - 1 (9)
где w0 - начальное условие задачи Коши, wi - решение системы ОДУ, определенное на ¿-м шаге, к - шаг, т - количество шагов. Координата t дискретизируется в сетку вида И = а + ¿к, где а - начальная координата рассматриваемого интервала.
Теперь рассмотрим неявный метод Эйлера, который, как и явный метод Эйлера, имеет 1 -й порядок точности. Запишем формулировку неявного метода Эйлера [1]:
w0 = а,
wi+1 = wi + hf(ti, wi+1), i = 0,1 ... m - 1 (10)
где w0 - начальное условие задачи Коши, wi - решение системы ОДУ, определенное на i-м шаге, h - шаг, т - количество шагов.
В данном случае свойство неявности вычислительного метода обусловлено наличием искомой величины wi+1 и в левой, и в правой частях уравнения (10).
Неизвестное значение вычисляется одним из методов решения нелинейных уравнений, в данном случае будет применяться метод root из модуля scipy.optimize. Аналогично предыдущему случаю, проверим функцию, реализующую неявный метод Эйлера. Результат представлен на рисунке 1.
I 40
- 30
г
- 20 10
Рисунок 1. Фазовый портрет траектории при г = 28 - аттрактор Лоренца (решение неявным методом Эйлера)
Теперь рассмотрим метод Адамса-Башфорта-Моултона. Это метод 4-го порядка точности. Формулировка метода Адамса-Башфорта-Моултона имеет
следующий вид [1]:
= а0, ш = а1, = а2, жЪ = а3, = а, (11)
Щ+1 [55/(t¿, ш;) - 59/(Ь-1, wi-1) + 37f{ti-2, wi-2) (12)
Щ+1 = [9f(ti+1, ^¿+1) + 19f(ti, wi) - 5/(Ь-1, wi-l) (13) + /(t¿-2,W¿-2)]
где I е 3 ... т.
Для определения начальных условий в данном случае использовалась функция 1трИсН_еи1ег.
Результат работы функции а(1ат5_Ьа5к[огМ_тоииоп представлен на рисунке 2.
Рисунок 2. Фазовый портрет траектории при г = 28 - аттрактор Лоренца (решение методом Адамса-Башфорта-Моултона)
ОПРЕДЕЛЕНИЕ ЧИСЛЕННОГО МЕТОДА ПО УМОЛЧАНИЮ
В качестве метода по умолчанию рекомендуется выбрать наиболее точный из трех методов. Очевидно, наибольшей точности можно добиться, используя метод Адамса-Башфорта-Моултона, т. к. он имеет 4-й порядок точности, в то время как явный и неявный методы Эйлера имеют 1 -й порядок точности.
Экспериментально было определено оптимальное значение временного шага - 0.001. Такой шаг позволяет получить довольно точное решение за адекватное время.
Также было экспериментально определено время - 100. Меньшее время позволило бы получить решение значительно быстрее, но для получения полного и точного фазового портрета лучше использовать время побольше.
АНАЛИЗ ДИНАМИЧЕСКИХ РЕЖИМОВ МОДЕЛИ ЛОРЕНЦА
Было сгенерировано 20 случайных начальных условий (этого достаточно для того, чтобы охарактеризовать любой из динамических режимов модели Лоренца, к тому же, выборка из небольшого количества начальных условий позволяет получить решение довольно быстро), исходя из равномерного распределения и следующих ограничений: х(0) е [-50; 50], у(0) е [-50; 50], г(0) е [0; 70].
Для каждого начального условия было получено решение ОДУ (в виде траектории) с использованием метода Адамса-Башфорта-Моултона.
Результат для г = 0 представлен на рисунке 3. Было принять решение строить фазовые портреты в 3d-формате, т. к. такие графики позволяют более
детально оценить поведение траектории.
Затем было изменено количество начальных условий (после всех итераций для всех г). Было сгенерировано всего 2 новых начальных условия -такие фазовые портреты позволяют более точно понять изменение траектории при различных г.
Рисунок 3. Фазовый портрет при г = 0 (для 20 начальных условий)
Можно сделать вывод, что при г = 0 присутствует сходимость решений в точке х = у = 2 = 0. Это не зависит от начальных условий и их количества. Аттрактором является начало координат, других устойчивых точек нет.
Рисунок 4. Фазовый портрет при г = 10 (для 20 начальных условий)
Траектории приходят в одну из двух стационарных точек: х = ^Ь(г — 1), у = ^Ь(г — 1), г = г — 1; х = —^Ь(г — 1), у = —^Ь(г — 1), г = г — 1
На рисунке 5 можно наблюдать ситуацию, когда две траектории пришли в одну стационарную точку (были сгенерированы новые случайные начальные условия).
Рисунок 6. Фазовый портрет при г = 30 (для 2 начальных условий)
Действительно, «странный аттрактор» появляется при г = 28.. .30.
Траектории будто бы поочередно спирально приближаются то к одной, то к другой стационарной точке. При этом возможны внезапные изменения траектории, которые не так легко предсказать. Именно при г = 28.30 система наиболее чувствительна к смене начальных условий: если во всех предыдущих случаях фазовые портреты при «старых» и «новых» начальных условиях были очень похожи друг на друга, то в данном случае фазовый портрет может значительно измениться при обновлении начальных условий.
Также рассмотрим фазовые портреты для 20 (рисунок 5) и 2 (рисунок 6) начальных условий при г = 100. В данном случае траектория претерпевает еще более серьезные изменения, чем в случае с г = 30. Скорее всего, при г = 100 (и более) система уже не подчиняется какому-либо математическому описанию, поэтому делать прогнозы о поведении траекторий практически невозможно. И, наверное, именно такая система лучше всего могла бы описать реальные механизмы, существующие в природе - случайные, зависящие от немыслимого количества различных факторов. Поэтому, несмотря на то, что решение при г = 100 (и более) довольно сложно понять и охарактеризовать, оно, пожалуй, является самым точным отражением реальных явлений. [4]
АНАЛИЗ ВЫБРАННОГО ДИНАМИЧЕСКОГО РЕЖИМА
Рассмотрим временные затраты, необходимые для нахождения решений при г = 28 (аттрактор Лоренца). Конечно, это не самый сложный динамический режим (гораздо более непредсказуемое и сложное поведение система демонстрирует при больших значениях г), но основной интерес в данной статье для нас представляет именно аттрактор Лоренца. Замерить время, необходимое для решения задачи, можно с помощью функции timeit. default timer(). Сначала определим время решения системы ОДУ для 20 начальных условий (временной интервал - 100, шаг по
умолчанию - 0.001, метод по умолчанию - метод Адамса-Башфорта-Моултона) - именно этот режим использовался для нахождения большей части траекторий (в предыдущем пункте). Будем увеличивать шаг по умолчанию в 2 раза (до 0.002) и уменьшать также в 2 раза (до 0.0005). Полученные результаты представлены в таблице 1.
Таблица 1. Временные затраты метода по умолчанию при различных Н (20 НУ)
Шаг (с) 0.002 0.001 0.0005
Затраченное время (с) 146.08 292.53 574.77
Результаты вполне соответствуют действительности: чем меньше шаг, тем больше вычислительных операций необходимо произвести, а значит, тем больше времени потребуется на решение задачи. [5]
Теперь сгенерируем 100 начальных условий. И определим время нахождения решений системы ОДУ для всех начальных условий различными методами: методом Эйлера, неявным методом Эйлера и методом Адамса-Башфорта-Моултона. В данном случае стоит уменьшить временной интервал (иначе потратится слишком много времени) - например, примем = 30. Будем рассматривать шаг по умолчанию 0.001, а также шаг 0.002 и шаг 0.0005. Временные затраты, необходимые для нахождения решений при использовании различных методов, представлены в таблице 2.
Таблица 2. Временные затраты рассматриваемых методов при различных Н (100 НУ)
Метод Метод Эйлера Неявный Метод
Шаг (с) метод Адамса-
Эйлера Башфорта-Моултона
0.002 22.23 с 304.98 с 227.21 с
0.001 38.77 с 552.53 с 424.87 с
0.0005 81.42 с 1077.91 с 826.42 с
Легко заметить, что самым быстрым методом является метод Эйлера. Это обусловлено относительной простотой операций, выполняемых на каждом шаге. Самым медленным методом является неявный метод Эйлера (в принципе, все неявные методы должны работать довольно медленно как раз в силу своей
«неявности» - требуется находить решение нелинейного уравнения, а эта операция порождает значительные временные затраты). Метод Адамса-Башфорта-Моултона тоже работает довольно медленно, так что высокой точности решения можно добиться только в ущерб скорости решения. Везде наблюдается увеличение времени нахождения решений при уменьшении шага, возникновение этого эффекта очевидно: чем меньше шаг, тем больше вычислительных операций требуется выполнить. Конечно, изменение шага также будет влиять на качество (точность и гладкость линий) фазового портрета.
Траектория, построенная методом Эйлера с шагом h = 0.01 является наименее точной, при детальном рассмотрении можно заметить, что это именно кусочно-линейная функция (видны переходы между отрезками). Траектории, построенные методом Эйлера с шагом h = 0.001 и шагом h = 0.0001 получились довольно точными. Визуально разница между этими решениями не так очевидна. Поэтому можно не дробить шаг далее - это увеличит время вычислений, но увидеть значительное улучшение точности фазового портрета не удастся.
Траектории, построенные с помощью метода Эйлера, неявного метода
Эйлера и метода Адамса-Башфорта-Моултона различаются, и это, скорее всего, обусловлено чувствительностью модели Лоренца к любым изменениям. Другие начальные условия, как и выбор численного метода (любой численный метод дает определенную погрешность) очень сильно влияют на вид фазового портрета.
ЗАКЛЮЧЕНИЕ
Были исследованы численные методы решения систем ОДУ (а именно: метод Эйлера, неявный метод Эйлера и метод Адамса-Башфорта-Моултона) на примере решения модели Лоренца.
Оказалось, что самым быстрым и простым методом решения систем ОДУ является метод Эйлера, однако этот метод значительно уступает в точности методу Адамса-Башфорта-Моултона. Неявный метод Эйлера требует наибольших временных затрат, однако полученное с использованием этого метода решение будет опять-таки менее точным, чем решение, полученное с помощью метода Адамса-Башфорта-Моултона. Пожалуй, неявный метод Эйлера не стоит использовать при решении подобных задач: если нужно получить решение быстро, то можно воспользоваться явным методом Эйлера, если в приоритете точность, то стоит использовать методы, например, 4-го порядка точности.
Также на качество полученного решения влияет выбор шага по времени: чем меньше шаг, тем точнее решение. Однако, визуальная оценка позволяет сделать предположение, что «бесконечно» дробить шаг не стоит, т. к. достаточно малое значение шага уже позволяет получить довольно точное решение (при дальнейшем уменьшении шага разница не очевидна). Но, чем меньше шаг, тем больше потребуется времени на решение, это связано с увеличением количества вычислительных операций. Поэтому выбор шага -это поиск компромисса между точностью и скоростью решения системы.
Список использованных источников
1. Першин А.Ю. Лекции по вычислительной математике. [Электронный ресурс] // Кафедра РК6 (Системы автоматизированного проектирования) МГТУ им. Н.Э. Баумана, 142 с. - (Дата обращения: 23.12.2020)
2. Динамические системы — сложное поведение. [Электронный ресурс] URL: http://w.ictnsc.ru/books/textbooks/akhmerov/ode/s-17/s- 17.html (Дата обращения: 23.12.2020)
3. Странные аттракторы [Электронный ресурс] URL: http://www.polybook.ru/comma/4.8.pdf (Дата обращения: 23.12.2020)
4. Модель Лоренца [Электронный ресурс] URL: https://keldysh.ru/pages/comma/html/ds/loren.htm (Дата обращения: 23.12.2020)
5. Ильяшенко Ю. С. Аттракторы и их фрактальная размерность. — М.: МЦНМО, 2005. — 16 с
Literature
1. Pershin A.Yu. Lectures on Computational Mathematics. [Electronic resource] // Department of RK6 (Computer-aided design systems) MGTU im. N.E. Bauman, 142 p. - (Date of access: 23.12.2020)
2. Dynamical systems - complex URL behavior: http://w.ict.nsc.ru/books/textbooks/akhmerov/ode/s-17/s- 17.html (Date of access: 2312.2020)
3. Strange attractors [Electronic resource] URL: http://www.polybook.ru/comma/4.8.pdf (Date of access: 23.12.2020)
4. Lorenz model [Electronic resource] URL: https://keldysh.ru/pages/comma/html/ds/loren.htm (Date of access: 23.12.202 0)
5. Ilyashenko Yu. S. Attractors and their fractal dimension. - M .: MTsNMO, 2005 .-- 16 p.
