Interested Article - Аттрактор Лоренца

решение системы при r =28 ― собственно, это и есть аттрактор Лоренца

решение системы при r =100 ― виден режим автоколебаний в системе

Аттрактор Лоренца (от англ. to attract — притягивать) ― странный аттрактор , впервые найденный Э. Н. Лоренцем в нелинейной системе обыкновенных дифференциальных уравнений

${\begin{cases}{\dot {x}}=\sigma (y-x)\\{\dot {y}}=x(r-z)-y\\{\dot {z}}=xy-bz\end{cases}}$

при следующих значениях параметров: σ=10, r =28, b =8/3. Эта система вначале была введена как первое нетривиальное галёркинское приближение для задачи о конвекции морской воды в плоском слое, чем и мотивировался выбор значений σ, r и b , но она возникает также и в других физических вопросах и моделях:

конвекция в замкнутой петле;
вращение водяного колеса;
модель одномодового лазера ;
диссипативный гармонический осциллятор с инерционной нелинейностью.

Исходная гидродинамическая система уравнений:

${\begin{cases}{\frac {\partial {\vec {v}}}{\partial t}}+\left({\vec {v}}\nabla \right){\vec {v}}=-{\frac {\nabla p}{\rho }}+\nu \nabla ^{2}{\vec {v}}+{\vec {g}}\\{\frac {\partial \rho }{\partial t}}+\nabla \cdot \left(\rho {\vec {v}}\right)=0\\{\frac {\partial T}{\partial t}}+\nabla \cdot \left(T{\vec {v}}\right)=\chi \nabla ^{2}T\\\rho =\rho _{0}\left(1-\gamma \left(T-T_{0}\right)\right)\end{cases}},$

где ${\vec {v}}$ — скорость течения, $T$ — температура жидкости, $T_{0}$ — температура верхней границы (на нижней поддерживается $T_{0}+\Delta T$ ), $\rho$ — плотность, $p$ — давление, ${\vec {g}}$ — сила тяжести, $\gamma ,\ \chi ,\ \nu$ — соответственно коэффициент теплового расширения , коэффициент температуропроводности и кинематической вязкости .

В задаче о конвекции модель возникает при разложении скорости течения и температуры в двумерные ряды Фурье и последующей их «обрезки» с точностью до первых-вторых гармоник. Кроме того, приведённая полная система уравнений гидродинамики записывается в приближении Буссинеска . Обрезка рядов в определённой мере оправдана, так как Сольцмен в своих работах продемонстрировал отсутствие каких-либо интересных особенностей в поведении большинства гармоник .

Применимость и соответствие реальности

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

Конвекция в плоском слое. Здесь x отвечает за скорость вращения водяных валов, y и z — за распределение температуры по горизонтали и вертикали, r — нормированное число Рэлея , σ — число Прандтля (отношение коэффициента кинематической вязкости к коэффициенту температуропроводности ), b содержит информацию о геометрии конвективной ячейки.
Конвекция в замкнутой петле. Здесь x — скорость течения, y — отклонение температуры от средней в точке, отстоящей от нижней точки петли на 90°, z — то же, но в нижней точке. Подведение тепла производится в нижней точке.
Вращение водяного колеса. Рассматривается задача о колесе, на ободе которого укреплены корзины с отверстиями в дне. Сверху на колесо симметрично относительно оси вращения льётся сплошной поток воды. Задача равнозначна предыдущей, перевернутой «вверх ногами», с заменой температуры на плотность распределения массы воды в корзинах по ободу.
Одномодовый лазер. Здесь x — амплитуда волн в резонаторе лазера, y — поляризация , z — инверсия населённостей энергетических уровней , b и σ — отношения коэффициентов релаксации инверсии и поля к коэффициенту релаксации поляризации, r — интенсивность накачки .

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

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

Примечательно везение Лоренца при выборе значения параметра $r$ , так как система приходит к странному аттрактору только при значениях, больших 24,74, при меньших поведение оказывается совершенно иным.

Поведение решения системы

Рассмотрим изменения в поведении решения системы Лоренца при различных значениях параметра r. На иллюстрациях к статье приведены результаты численного моделирования для точек с начальными координатами (10,10,10) и (-10,-10,10). Моделирование производилось с помощью приведённой ниже программы, написанной на языке Фортран , построение графиков по полученным таблицам — из-за слабых графических возможностей Фортрана с помощью Compaq Array Viewer.

r <1 — аттрактором является начало координат, других устойчивых точек нет.

1< r <13,927 — траектории спирально приближаются (это соответствует наличию затухающих колебаний) к двум точкам, положение которых определяется формулами:

${\begin{cases}x=\pm {\sqrt {b(r-1)}}\\y=\pm {\sqrt {b(r-1)}}\\z=r-1\end{cases}}$

Эти точки определяют состояния стационарного режима конвекции, когда в слое формируется структура из вращающихся валов жидкости.

r ≈13,927 — если траектория выходит из начала координат, то, совершив полный оборот вокруг одной из устойчивых точек, она вернется обратно в начальную точку — возникают две гомоклинические петли. Понятие гомоклинической траектории означает, что она выходит и приходит в одно и то же положение равновесия.

r >13,927 — в зависимости от направления траектория приходит в одну из двух устойчивых точек. Гомоклинические петли перерождаются в неустойчивые предельные циклы, также возникает семейство сложно устроенных траекторий, не являющееся аттрактором, а скорее наоборот, отталкивающее от себя траектории. Иногда по аналогии эта структура называется «странным репеллером» (англ. to repel — отталкивать).

r ≈24,06 — траектории теперь ведут не к устойчивым точкам, а асимптотически приближаются к неустойчивым предельным циклам — возникает собственно аттрактор Лоренца. Однако обе устойчивые точки сохраняются вплоть до значений r ≈24,74.

r ≈28 — классическое значение параметра, рассмотренное в статье Лоренца. Все три состояния равновесия являются неустойчивыми и траектории из их окрестностей притягиваются к хаотическому (локальному) аттрактору (который, таким образом, является относительно всех состояний равновесия). Хаотический аттрактор имеет дробную , для которой аналитическая оценка сверху может быть получена аналитически через форму ляпуновской размерности глобального аттрактора, а оценка снизу может быть получена аналитико-численно через ляпуновскую размерность неустойчивых периодических траекторий на аттракторе . Приближения к таким траекториям могут быть найдены с высокой точностью методом гармонического баланса . Для высокоточного численного моделирования динамики системы Лоренца обычно используется метод степенных рядов .

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

Значимость модели

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

Из-за своей характерной формы аттрактор получил название «бабочка Лоренца», что породило понятие « эффект бабочки » в теории хаоса , впоследствии ошибочно связанное в массовом сознании с известным рассказом Рэя Брэдбери .

Программы, моделирующие поведение системы Лоренца

Borland C

#include <graphics.h>
#include <conio.h>
void main()
{
    double x = 3.051522, y = 1.582542, z = 15.62388, x1, y1, z1;
    double dt = 0.0001;
    int a = 5, b = 15, c = 1;
    int gd=DETECT, gm;
    initgraph(&gd, &gm, "C:\\BORLANDC\\BGI");
    do {
	x1 = x + a*(-x+y)*dt;
	y1 = y + (b*x-y-z*x)*dt;
	z1 = z + (-c*z+x*y)*dt;
	x = x1;	y = y1;	z = z1;
	putpixel((int)(19.3*(y - x*0.292893) + 320),
		 (int)(-11*(z + x*0.292893) + 392), 9);
    } while (!kbhit());
    closegraph();
}

Mathematica

data = Table[
   With[{N = 1000, dt = 0.01, a = 5, b = 1 + j, c = 1},
    NestList[Module[{x, y, z, x1, y1, z1},
       {x, y, z} = #;
       x1 = x + a (-x + y) dt;
       y1 = y + (b x - y - z x) dt;
       z1 = z + (-c z + x y) dt;
       {x1, y1, z1}] &,
     {3.051522, 1.582542, 15.62388}, N
     ]
    ],
   {j, 0, 5}];
Graphics3D@MapIndexed[{Hue[0.1 First[#2]], Point[#1]} &, data]

JavaScript и HTML5

<html>
<body>
  <canvas height='500' width='500' id='cnv'></canvas>
  <script>
        var cnv = document.getElementById("cnv");
        var cx = cnv.getContext('2d');
        var x = 3.051522, y = 1.582542, z = 15.62388, x1, y1, z1;
        var dt = 0.0001;
        var a = 5, b = 15, c = 1;
        var h = parseInt(cnv.getAttribute("height"));
        var w = parseInt(cnv.getAttribute("width"));
        var id = cx.createImageData(w, h);
        var rd = Math.round;
        var idx = 0;
        i=1000000; while (i--) {
                x1 = x + a*(-x+y)*dt;
                y1 = y + (b*x-y-z*x)*dt;
                z1 = z + (-c*z+x*y)*dt;
                x = x1; y = y1; z = z1;                         
                idx=4*(rd(19.3*(y - x*0.292893) + 320) + rd(-11*(z + x*0.292893) + 392)*w);
                id.data[idx+3] = 255;
        }
        cx.putImageData(id, 0, 0);
  </script>
</body>
</html>

MATLAB

%Solution for the Lorenz equations in the time interval [0,100] with initial conditions [1,1,1].
clear all
clc
sigma=10;
beta=8/3;
rho=28;
f = @(t,a) [-sigma*a(1) + sigma*a(2); rho*a(1) - a(2) - a(1)*a(3); -beta*a(3) + a(1)*a(2)];
%'f' is the set of differential equations and 'a' is an array containing values of x,y, and z variables.
%'t' is the time variable
[t,a] = ode45(f,[0 100],[1 1 1]);%'ode45' uses adaptive Runge-Kutta method of 4th and 5th order to solve differential equations
plot3(a(:,1),a(:,2),a(:,3)) %'plot3' is the command to make 3D plot

Maxima

 --> load(dynamics)$
 [sigma, r,b]: [10,28,8/3]$
 eq: [sigma*(y-x), x*(r-z)-y, x*y-b*z]$
 init: [1.0,0,0]$
 t_range: [t,0,50,0.01]$
 sol: rk(eq, [x, y,z], init, t_range)$
 len: length(sol)$
 t: makelist(sol[k][1], k,1,len)$
 x: makelist(sol[k][2], k,1,len)$
 y: makelist(sol[k][3], k,1,len)$
 z: makelist(sol[k][4], k,1,len)$
 plot2d([discrete, t,x])$
 
 --> load(draw)$
 draw3d(point_size=0.01,  points_joined=true,
   point_type=filled_circle,points(x,y,z)
    )$

Python

"""
================
Lorenz Attractor
================
"""
import numpy as np
import matplotlib.pyplot as plt

def lorenz(x, y, z, s=10, r=28, b=2.667):
    '''
    Given:
       x, y, z: a point of interest in three dimensional space
       s, r, b: parameters defining the lorenz attractor
    Returns:
       x_dot, y_dot, z_dot: values of the lorenz attractor's partial
           derivatives at the point x, y, z
    '''
    x_dot = s*(y - x)
    y_dot = r*x - y - x*z
    z_dot = x*y - b*z
    return x_dot, y_dot, z_dot

dt = 0.01
num_steps = 10000

# Need one more for the initial values
xs = np.empty(num_steps + 1)
ys = np.empty(num_steps + 1)
zs = np.empty(num_steps + 1)

# Set initial values
xs[0], ys[0], zs[0] = (0., 1., 1.05)

# Step through "time", calculating the partial derivatives at the current point
# and using them to estimate the next point
for i in range(num_steps):
    x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])
    xs[i + 1] = xs[i] + (x_dot * dt)
    ys[i + 1] = ys[i] + (y_dot * dt)
    zs[i + 1] = zs[i] + (z_dot * dt)


# Plot
fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot(xs, ys, zs, lw=0.5)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")

plt.savefig('Lorenz Attractor')
plt.show()

Примечания

Saltzman, Barry (1962). «Finite Amplitude Free Convection as an Initial Value Problem—I». Journal of the Atmospheric Sciences 19 (4): 329—341.
Kuznetsov, N.V.; Mokaev, T.N.; Kuznetsova, O.A.; Kudryashova, E.V. (2020). . Nonlinear Dynamics . doi : . из оригинала 28 июня 2021 . Дата обращения: 20 сентября 2020 .
Leonov, G.A.; Kuznetsov, N.V.; Korzhemanova, N.A.; Kusakin, D.V. (2016). "Lyapunov dimension formula for the global attractor of the Lorenz system". Communications in Nonlinear Science and Numerical Simulation . 41 : 84—103. arXiv : . Bibcode : . doi : . {{ cite journal }} : Недопустимый |ref=harv ( справка )
Kuznetsov, Nikolay. / Nikolay Kuznetsov, Volker Reitmann. — Cham : Springer, 2021. . Дата обращения: 20 сентября 2020. Архивировано 3 июня 2020 года.
Pchelintsev, A.N. (2020). . Дифференциальные уравнения и процессы управления (4): 59—75. arXiv : .
Pchelintsev, A.N. (2014). "Numerical and physical modeling of the dynamics of the Lorenz system". Numerical Analysis and Applications . 7 (2): 159—167. doi : . S2CID .

Литература

Кузнецов С. П. , Лекция 3. Система Лоренца; Лекция 4. Динамика системы Лоренца. // — М.: Физматлит, 2001.
Saltzman B . Finite amplitude free convection as an initial value problem. // Journal of the atmospheric science, № 7, 1962 — p. 329—341.
Лоренц Э . Детерминированное непериодическое движение // Странные аттракторы. — М., 1981. — С. 88-116.