Аттрактор Лоренца

решение системы при r=0,3

решение системы при r=1,8

решение системы при r=3,7

решение системы при r=10

решение системы при r=16

решение системы при r=24,06

решение системы при 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.add_subplot(111, 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()

Литература

• Кузнецов С. П., Лекция 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.

См. также

• Динамический хаос
• Теория хаоса
• Лемниската Бернулли
• Васцилляция