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

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

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

Дивний атрактор Лоренца при $\lambda =28,\sigma =10,b=8/3$.

Рівняння

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

${\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}},$

${\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}}~,$

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

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

${\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),$
${\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),$

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

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

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

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

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

Стаціонарні точки

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

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

Траєкторії

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

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

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

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

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

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

${\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).$

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

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

Програми, які моделюють поведінку системи рівнянь Лоренца

Matlab

% 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))

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

Borland Pascal

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.

Python

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()

FORTRAN

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

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]

QBASIC/FreeBASIC(«fbc -lang qb»)

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.]. (англ.)