Дивний атрактор Лоренца

Ди́вний атра́ктор Ло́ренца — атрактор, що демонструє хаотичну поведінку і є розв'язком системи трьох нелінійних диференціальних рівнянь, вперше записаних^[1] в 1963 році Едвардом Лоренцом при розгляді конвекційного руху в однорідному шарі рідини, що підігрівається знизу. Рівняння Лоренца також описують конвекцію в кільцевій трубці^[2] та поведінку одномодового лазера. Належить до класу так званих дивних атракторів. Варто зазначити, терміни хаос та дивний атрактор не вживалися в оригінальній роботі Лоренца (вони з'явилися в науковій літературі дещо пізніше), натомість йшлося про аперіодичні рухи.

Початковою системою, яка в кінцевому результаті призводить до атрактора Лоренца, є однорідний шар рідини висотою H та з фіксованою різницею температур, ΔT між верхнім та нижнім рівнями. Якщо припустити, що всі рухи рідини паралельні площині XZ, та не відбувається змін в напрямку осі Y, то, записавши рівняння Нав'є — Стокса, рівняння неперервності рідини, рівняння теплопровідності та скориставшись наближенням Бусінеска, можна отримати рівняння руху рідини у такій формі:

${\displaystyle {\frac {\partial }{\partial t}}\nabla ^{2}\psi ={\frac {\partial \psi }{\partial z}}{\frac {\partial \nabla ^{2}\psi }{\partial x}}-{\frac {\partial \psi }{\partial x}}{\frac {\partial \nabla ^{2}\psi }{\partial z}}+\nu \nabla ^{2}(\nabla ^{2}\psi )+g\alpha {\frac {\partial T}{\partial x}},}$

${\displaystyle {\frac {\partial T}{\partial t}}={\frac {\partial T}{\partial z}}{\frac {\partial \psi }{\partial x}}-{\frac {\partial T}{\partial x}}{\frac {\partial \psi }{\partial z}}+\kappa \nabla ^{2}T+{\frac {\Delta T}{H}}{\frac {\partial \psi }{\partial x}}~,}$

де ${\displaystyle \psi }$ — функція потоку рідини (${\displaystyle v_{x}=-\partial _{z}\psi }$, ${\displaystyle v_{z}=\partial _{x}\psi }$, ${\displaystyle {\vec {v}}=(v_{x},v_{z})}$ — поле швидкостей), ${\displaystyle T}$ — відхилення температури рідини від її критичного значення, за якого зникає конвекція. Параметри ${\displaystyle g,\alpha ,\nu ,\kappa }$ позначають, відповідно, прискорення вільного падіння, коефіцієнт теплового розширення, кінематичну в'язкість рідини та її теплопровідність. Як було встановлено Релеєм, у такій системі можуть утворюватися конвекційні вали (ізотропні в напрямку осі OY), у яких відбувається коловорот рідини: тепліша рідина піднімається нагору, а холодніша опускається донизу. Конвекційні вали ділять площину XZ на приблизно однакові комірки.

Згідно з ідеями Релея^[3] та Зальцмана^[4], можна розкласти поля ${\displaystyle \psi (x,z,t)}$ та ${\displaystyle T(x,z,t)}$ в ряд Фур'є по ${\displaystyle x}$ та ${\displaystyle z}$ і обмежитись лише першими членами розкладу:

${\displaystyle {\frac {a}{\kappa (1+a^{2})}}\psi (x,z,t)={\sqrt {2}}X\sin \left({\frac {\pi ax}{H}}\right)\sin \left({\frac {\pi z}{H}}\right),}$

${\displaystyle {\frac {\pi R_{a}}{\Delta TR_{c}}}T(x,z,t)={\sqrt {2}}Y\cos \left({\frac {\pi ax}{H}}\right)\sin \left({\frac {\pi z}{H}}\right)-Z\sin \left({\frac {2\pi z}{H}}\right),}$

де ${\displaystyle R_{a}=g\alpha H^{3}\Delta T/(\nu \kappa )}$ — число Релея, ${\displaystyle R_{c}=\pi ^{4}(1+a^{2})^{3}/a^{2}}$, параметр ${\displaystyle a}$ задає співвідношення вертикального та горизонтального розмірів комірок, що утворюються в результаті конквекції, а змінні ${\displaystyle X,Y,Z}$ залежать лише від часу. Після переходу до безрозмірного часу шляхом заміни ${\displaystyle t\to t\cdot \pi ^{2}(1+a^{2})\kappa /H^{2}}$ отримується система трьох звичайних диференціальних рівнянь, що носить назву рівнянь Лоренца:

${\displaystyle {\dot {X}}=-\sigma X+\sigma Y,}$

${\displaystyle {\dot {Y}}=-XZ+\lambda X-Y,}$

${\displaystyle {\dot {Z}}=XY-bZ,}$

де точка означає диференціювання за часом, ${\displaystyle \sigma =\nu /\kappa }$ — число Прандтля, ${\displaystyle \lambda =R_{a}/R_{c}}$ та ${\displaystyle b=4/(1+a^{2})}$. Динамічні змінні ${\displaystyle X(t)}$, ${\displaystyle Y(t)}$ та ${\displaystyle Z(t)}$ описують відповідно інтенсивність конвективного руху, різницю температур висхідного та низхідного потоків рідини та відхилення вертикального розподілу температури від лінійного режиму.

Аналіз властивостей стаціонарних точок системи зручно робити, змінюючи параметр ${\displaystyle \lambda }$:

• ${\displaystyle 0<\lambda <1}$: стаціонарна точка ${\displaystyle (0,0,0)}$ відповідає стану відсутності конвекції, при ${\displaystyle 0<\lambda <1}$ вона є стійкою (стійкий вузол), при ${\displaystyle \lambda >1}$ стає нестійкою (сідло-вузлом);
• ${\displaystyle 1<\lambda <\lambda _{1}}$: в момент втрати стійкості точки ${\displaystyle (0,0,0)}$ з'являються дві інші стійкі стаціонарні точки — фокуси: ${\displaystyle ({\sqrt {b(\lambda -1)}},{\sqrt {b(\lambda -1)}},\lambda -1)}$ та ${\displaystyle (-{\sqrt {b(\lambda -1)}},-{\sqrt {b(\lambda -1)}},\lambda -1)}$, які відповідають режиму сталої конвекції. Фазові криві прямують до цих стаціонарних точок по спіралях. Чим більшим є параметр ${\displaystyle \lambda }$, тим більший розмах мають ці спіралі при підході до стаціонарних точок;
• ${\displaystyle \lambda _{1}<\lambda <\lambda _{2}}$: це критичне значення ${\displaystyle \lambda _{1}}$ можна встановити лише чисельно, зокрема для ${\displaystyle \sigma =10,b=8/3}$ воно дорівнює ${\displaystyle \lambda _{1}\simeq 13,926}$. При ${\displaystyle \lambda =\lambda _{1}}$ траєкторія, почавши рух з початку координат ${\displaystyle (0,0,0)}$ знову приходить в нього. Таким чином, відбувається нелокальна біфуркація. При ${\displaystyle \lambda >\lambda _{1}}$ навколо стаціонарних точок з'являються два нестійких граничних цикли. Утворюється інваріантна множина точок, що відповідає хаотичному «блуканню» траєкторій, що постійно відштовхуються від граничних циклів. Проте ця множина не є притягуючою, тому говорять про дивний репелер;
• ${\displaystyle \lambda _{2}<\lambda <{\frac {\sigma (\sigma +b+3)}{\sigma -b-1}}}$: значення ${\displaystyle \lambda _{2}}$ (при ${\displaystyle \sigma =10,b=8/3}$ воно дорівнює ${\displaystyle \lambda _{2}\simeq 24,06}$) сигналізує про перетворення дивного репелера в дивний атрактор (іноді говорять про т. зв. нестандартний атрактор Лоренца), який співіснує з двома іншими атракторами — стійкими фокусами;
• ${\displaystyle \lambda >{\frac {\sigma (\sigma +b+3)}{\sigma -b-1}}}$: при досягненні цього критичного значення нестійкі граничні цикли стягуються в стаціонарні точки, і стаціонарні точки втрачають свою стійкість; утворюється т. зв. стандартний атрактор Лоренца, що стає єдиним атрактором системи.

Фазові траєкторії, що починаються в будь-якій точці, спочатку, здавалося б, притягаються до однієї зі стаціонарних точок, але не можуть підійти надто близько, оскільки стаціонарна точка нестійка. У якийсь момент фазова траєкторія перестрибує в окіл іншої стаціонарної точки, але там теж не може затриматися, й знову перестрибує до першої стаціонарної точки. Як наслідок, система безупинно аперіодично перестрибує від однієї точки до іншої.

Фазові траєкторії чутливі до щонайменшої зміни початкових умов. Дві нескінченно близькі початкові точки в фазовому просторі з часом розходяться.

Ефект метелика
Момент часу t = 1	Момент часу t = 2	Момент часу t = 3
Рисунки зображають поведінку двох траєкторій рівнянь Лоренцца (при параметрах ${\displaystyle \lambda =28,\sigma =10,b=8/3}$), які в початковий момент були відокремлені одна від одної по змінній ${\displaystyle X}$ на величину 10−5. В початковий момент часу складається враження, що траєкторії збігаються, проте через деякий час стає очевидним, що синя та жовта траєкторії суттєво розбігаються.
Java анімація атрактора Лоренца. [Архівовано 11 березня 2008 у Portugese Web Archive]

У 1983 році Грассбергер та Прокаччія оцінили^[5] розмірність Гаусдорфа дивного атрактора Лоренца й одержали величину 2.06 ± 0.01.

При дуже великих значеннях ${\displaystyle \lambda }$, ${\displaystyle \lambda \gg \sigma (\sigma +b+3)/(\sigma -b-1)}$ динаміка системи Лоренца описується звичайними граничними циклами. Зокрема, при ${\displaystyle \lambda =99,96,\sigma =10,b=8/3}$ утворюється граничний цикл, що має вигляд вузлового тора ${\displaystyle T(3,2)}$. При зменшенні ${\displaystyle \lambda }$ перехід до хаотичного режиму відбувається чарез каскад біфуркацій подвоєння періоду.

Водночас фазові траєкторії не можуть втекти на нескінченність, оскільки при великих X, Y, Z виникають сили, що повертають фазові траєкторії в область малих значень змінних. Якщо домножити перше рівняння Лоренца на ${\displaystyle X/\sigma }$, друге — на ${\displaystyle Y}$, а третє — на ${\displaystyle Z}$, а потім додати всі три рівняння, то результат можна записати як:

${\displaystyle {\frac {1}{2}}{\frac {d}{dt}}\left({\frac {X^{2}}{\sigma }}+Y^{2}+Z^{2}\right)=-\left(X-{\frac {Y}{2}}\right)^{2}-{\frac {3}{4}}Y^{2}-b\left(Z-{\frac {\lambda }{2}}\right)^{2}+{\frac {b\lambda ^{2}}{4}}\equiv E(X,Y,Z).}$

Поверхня, що задається нерівністю ${\displaystyle E(X,Y,Z)\geq 0}$, має вигляд еліпсоїда зі зміщеним центром мас. Нескладно здогадатися, що за будь-якого вибору початкових умов еволюція системи на атракторі не призведе до виходу за межі одного з таких еліпсоїдів.

Квазіперіодичних коливань в системі Лоренца не може бути за жодних умов.

% Solve over time interval [0,100] with initial conditions [1,1,1]
% ''f'' is set of differential equations
% ''a'' is array containing x, y, and z variables
% ''t'' is time variable

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)];
[t,a] = ode45(f,[0 100],[1 1 1]);     % Runge-Kutta 4th/5th order ODE solver
plot3(a(:,1),a(:,2),a(:,3))

#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();
}

Program Lorenz;
Uses CRT, Graph;
Const
  x: Real = 3.051522;
  y: Real = 1.582542;
  z: Real = 15.62388;
  dt = 0.0001;
  a = 5;
  b = 15;
  c = 1;
Var
  gd, gm: Integer;
  x1, y1, z1: Real;
Begin
  gd:=Detect;
  InitGraph(gd, gm, 'c:\bp\bgi');
  While not KeyPressed Do Begin
      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(Round(19.3*(y - x*0.292893) + 320),
               Round(-11*(z + x*0.292893) + 392), 9);
    End;
    CloseGraph;
    ReadKey;
End.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D

rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0

def f(state, t):
    x, y, z = state  # Unpack the state vector
    return sigma * (y - x), x * (rho - z) - y, x * y - beta * z  # Derivatives

state0 = [1.0, 1.0, 1.0]
t = np.arange(0.0, 40.0, 0.01)

states = odeint(f, state0, t)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(states[:, 0], states[:, 1], states[:, 2])
plt.show()

program LorenzSystem

real,parameter::sigma=10
real,parameter::r=28
real,parameter::b=2.666666
real,parameter::dt=.01
integer,parameter::n=1000

real x,y,z

open(1,file='result.txt',form='formatted',status='replace',action='write')

x=10.;y=10.;z=10.

do i=1,n,1
    x1=x+sigma*(y-x)*dt
    y1=y+(r*x-x*z-y)*dt
    z1=z+(x*y-b*z)*dt
    x=x1
    y=y1
    z=z1
    write(1,*)x,y,z
enddo

print *,'Done'

close(1)

end program LorenzSystem

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]

DIM x, y, z, dt, x1, y1, z1 AS SINGLE
DIM a, b, c AS INTEGER
x = 3.051522: y = 1.582542: z = 15.62388: dt = 0.0001
a = 5: b = 15: c = 1
SCREEN 12
PRINT "Press Esc to quit"
WHILE INKEY$ <> CHR$(27)
    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
    PSET ((19.3 * (y - x * .292893) + 300), (-11 * (z + x * .292893) + 360)), 9
WEND
END

Вікісховище має мультимедійні дані за темою: Дивний атрактор Лоренца

• Сугаков В. Й. Основи синерґетики. — К. : Обереги, 2001. — 287 с.
• Кузнецов С. П. Динамический хаос. — М. : Физматлит, 2006. — 256 с.
• Лихтенберг А., Либерман М. Регулярная и стохастическя динамика. — М. : Мир, 1984. — 528 с.
• Лоренц Э. Детерминированное непериодическое течение //  — М. : Мир, 1981. — С. 88—116.
• Шустер Г. Детерминированный хаос. Введение. — М. : Мир, 1988. — 248 с.
• Ott E. Chaos in Dynamical Systems. — Cambridge University Press, 2002. — ISBN 0-521-01084-5.
• Weisstein E. W. Lorenz attractor [Архівовано 9 квітня 2008 у Wayback Machine.]. (англ.)