1. Причини інтересу до системи Лоренца
У даній публікації розглянуто складну поведінку динамічної системи на прикладі потоку, тобто зсуву за часом t вздовж траекторій — розв'язків системи звичайних диференціальних рівнянь:
| dx1/dt = – σx1 + σx2 ; | (1) |
| dx2/dt = – x1x3 + rx1 – x2 ; | (2) |
| dx3/dt = x1x2 – bx3 . | (3) |
Тут σ, b та r — деякі додатні параметри.
Cистему (1–3) запровадив Е.Лоренц у ході моделювання конвектиного потоку рідини, що підігрівають знизу (Edward Norton Lorenz. Deterministic Nonperiodic Flow // Journal of the Atmospheric Sciences, 1963, volume 20, issue 2 (March), pp. 130–141). Такий потік описують за допомогою системи диференціальних рівнянь у часткових похідних. Систему Лоренца (1–3) отримують за допомогою певного наближення — проектування на деякий тривимірний підпростір нескінченовимірного простору функцій. Інакше кажучи, використанням методу Бубнова (1913) — Гальоркіна (Галёркин Б. Г. Стержни и пластинки. Ряды в некоторых вопросах упругого равновесия стержней и пластинок. // Вестник инженеров. — 1915. — Т. 1. — С. 897—908).
У ході чисельного інтегрування системи (1–3) Е.Лоренц виявив, що при σ = 10, b = 8/3, r = 28 траекторії цієї динамічної системи одночасно демонструють хаотичну нерегулярну поведінку і притягуються до певної множини зі складною структурою, яку називають аттрактором (від. англ. to attract — притягувати). На рис. 1–3 подано залежність від часу t∈[0; 20] розв'язків системи рівнянь при таких початкових величин:
при вказаних величинах параметрів сисеми Лоренца. Ці величини отримано інтегруванням системи Лоренца при тих самих величинах параметрів на проміжку [0; 106] і таких початкових величинах: x1 = x2 = x3 = 1.
Такою поведінку траекторій природно пов’язали з турбулентним потоком рідини. Але остаточно невідомо, чи турбулентність пов'язана з хаотичністю руху у скінченовимірному просторі коефіцієнтів проекції, чи саме з нескінченовимірністю простору функцій, що описують потік рідини.
Звичайно, пересічні люди далекі від проблем опису турбулентності. Але прогнозами погоди цікавляться усі, у тому числі, з економічних міркуваню. Система Лоренца виникла саме у процесі математичного моделювання метеорологічних процесів з метою створення точних синоптичних прогнозів. Сучасні моделі метеорологічних процесів ґрунтуються на нелінійних диференціальних рівняннях у часткових похідних. Можливо, непередбачуваність погоди пояснюється наявністю складної структури відповідного аттрактора, а не лише випадковими втручанням космічних сил.
Питання має значну передісторію. Ще у 1814 році французький вчений П’єр-Симон Лаплас запровадив поняття демона, що був предметом наукових дискусій протягом тривалого часу. Видуманий демон знав розташування і швидкості кожної частинки у Всесвіті і, знаючи усі фізичні закони, міг передбачити мабутнє кожної частинки й описати її минуле. Чи можливий такий демон теоретично? Успіхи науки Нового часу (розрахунок орбіт планет, передбачення появи комет) давали надію на ствердну відповідь.
У 1927 році німецький фізик Вернер Гайзенберг відкрив принцип невизначенності. Пояснимо принцип на прикладі пари декартова координата + проекція імпульса на відповідний напрям. Розглянемо ансамбль (сукупність) еквівалентних частинок, що не взаємодіють між собою і знаходяться в одному стані. При вимірюванні координати x та імпульса p результати вимірювання будуть випадковими, а їхні дисперсії (середньо квадратичні відхилення від середніх значень) Δx та Δp задовольняють таку нерівність:
Δx · Δp ≥ ћ /2.
Тут ћ = 1,054 571 628(53) · 10–34 Дж·с — стала Дірака (зведена стала Планка). У загальному випадку принцип невизначеності виникає між довільними змінними, яким відповідають оператори, що не комутують (не переставляються). Таким чином, квантові системи принципово непідвладні демону Лапласа.
Згодом зауважили, що існування демона суперечило б законам термодинамики: йому не вистачило б знань та інформаційних ресурсів Всесвіту для вичерпного передбачення його майбутнього. Але чи є принципові труднощі для демона Лапласа (читай суперкоп’ютера) у випадку ізольованої скінченої повністю детермінованої системи?
Довільні фізичні вимірювання містять певну помилку. Змінні у пам’яті комп'ютера зберігаються з певною точністю.
Можливо достатньо обмежитися кількома значащими цифрами у початкових данних, щоб бути впевненим у кінцевому результаті хоча б з певною точністю?
Система Лоренца — це хронологічно перший приклад системи звичайних диференціальних рівнянь, розв’язки якої суттєво чутливі до початкових значень. Настільки чутливі, що породили поняття «детермінований хаос», яке традиційно вживають для позначення нерегулярної та непередбачуваної поведінки детермінованих нелінійних динамічних систем. Чи округлення у теоретичній моделі, чи похибка вимірювання у реальній системі — і система веде себе кардинально по іншому.
Лоренц міркував так: якщо погода справді належить до настільки чутливих систем, то помах крил чайки може викликати суттєві зміни погоди. Згодом чайку замінили на метелика, а в 1972 році з’явилася работа «Передбачуваність: чи може помах крил метелика в Бразилії призвести до торнадо в Техасі?». Так народився знаменитий термін «ефект метелика», що викликає асоціації з науково-фантастичним оповіданням Рея Бредбері та з відкриттям Лоренца – дивним аттрактором, названим на його честь.
2. Результати часткового дослідження системи Лоренца
Стаціонарні (нерухомі) точки системи Лоренца можна знайти, прирівнявши до нуля праві частини рівнянь (1–3):
Матриця ∂ f (x) / ∂ x (похідна подя напрямів) для системи Лоренца має такий вигляд.
| – σ | σ | 0 |
| r – x3 | – 1 | – x1 |
| x2 | x1 | – b |
Власні числа λ матриці ∂ f (x) /∂ x задовольняють характеристичне рівняння:
det (λ – ∂ f (x) /∂ x) = 0.
Тут det (A) — визначник матриці A.
Подамо сценарій зміни поведінки системи Лоренца (1–3) при сталих σ = 10, b = 8/3 і зростанні r від нуля.
Аналітичні дослідження стійкості нерухомих точок можна проводити і традиційно «на кінчику пера», і за допомогою спеціального програмного забезпечення. Скористаємося нагодою, щоб показати використання системи алгебричних обчислень Reduce, яку можна разом з описом запозичити з такого сайту: http://reduce-algebra.sourceforge.net/. Код програми має такий вигляд (знак % позначає початок коментаря, кінець якого є кінцем відповідного рядка).
off nat, % Виведення в один рядок
echo; % Блокування виведення ракції на команди
order z; % Запис характеристичного рівняння
% відносно власного числа z
% у порядку спадання степенів z
factor z; % Групування одночленів з однаковим
% степенем z
n:=3; % Розмірність фазового простору -
% кількість динамічних змінних
% Виділення назв для:
array x(n), % динамічних змінних х,
f(n), % поля напрямків f,
l(n); % власних чисел fx
matrix fx(n,n); % матриці похідних f за x
% Задання вектор-стовпчиха динамічних змінних
x(1):=x1;
x(2):=x2;
x(3):=x3;
% Задання поля напрямків f
f(1):=s*(x2-x1);
f(2):=-x1*x3+r*x1-x2;
f(3):= x1*x2-b*x3;
% Визначення значень параметрів
s:=10;
b:=8/3;
% Знаходження похідної поля напрямків
for j:=1:n do
for k:=1:n do fx(j,k):=df(f(j),x(k));
for j:=1:n do fx(j,j):=fx(j,j)-z;
% Знаходження списку стаціонарних точок x0
x0:=solve(
{f(1),f(2),f(3)}, % Ліві частини рівнянь за умови,
% що праві частини дорівнюють 0
{x1,x2,x3}); % невідомі змінні
n0:=length(x0); % Кількість стаціонарних точок
% (у тому числі уявних)
out "result.txt"; % Визначення файлу виведення
for j:=1:n0 do <<
write("______________________________________");
write("Координати стаціонарної точки № ",j,":");
for k:=1:n do <<
y:=part(x0,j,k);
write(y);
let y >>;
write("Похідна поля напрямків:");
let z=0;
for k:=1:n do for l:=1:n do
write("df",k," / dx",l," = ",fx(k,l));
clear z;
polynom:=det(fx);
write("Характеристичне рівняння: ",
(-1)**n*polynom," = 0");
z0:=solve(polynom=0,z); % Cписок власних чисел
if length(z0)>1 then <<
write("Власні значення - корені ",
"характеристичного рівняння:");
for k:=1:n do write(part(z0,k)) >>;
clear x1,x2,x3 % Вивільнення змiних;
>>;
shut "result.txt";% Закриття файлу виведення
end;Після виконання цієї програми файл result.txt матиме такий вигдяд (вилучено порожні рядки).
______________________________________ Координати стаціонарної точки № 1: x1=(2*sqrt(r - 1)*sqrt(2))/sqrt(3) x2=(2*sqrt(r - 1)*sqrt(2))/sqrt(3) x3=r - 1 Похідна поля напрямків: df1 / dx1 = -10 df1 / dx2 = 10 df1 / dx3 = 0 df2 / dx1 = 1 df2 / dx2 = -1 df2 / dx3 = ( - 2*sqrt(r - 1)*sqrt(2))/sqrt(3) df3 / dx1 = (2*sqrt(r - 1)*sqrt(2))/sqrt(3) df3 / dx2 = (2*sqrt(r - 1)*sqrt(2))/sqrt(3) df3 / dx3 = ( - 8)/3 Характеристичне рівняння: (3*z**3 + 41*z**2 + 8*z*(r + 10) + 160*(r - 1))/3 = 0 ______________________________________ Координати стаціонарної точки № 2: x1=( - 2*sqrt(r - 1)*sqrt(2))/sqrt(3) x2=( - 2*sqrt(r - 1)*sqrt(2))/sqrt(3) x3=r - 1 Похідна поля напрямків: df1 / dx1 = -10 df1 / dx2 = 10 df1 / dx3 = 0 df2 / dx1 = 1 df2 / dx2 = -1 df2 / dx3 = (2*sqrt(r - 1)*sqrt(2))/sqrt(3) df3 / dx1 = ( - 2*sqrt(r - 1)*sqrt(2))/sqrt(3) df3 / dx2 = ( - 2*sqrt(r - 1)*sqrt(2))/sqrt(3) df3 / dx3 = ( - 8)/3 Характеристичне рівняння: (3*z**3 + 41*z**2 + 8*z*(r + 10) + 160*(r - 1))/3 = 0 ______________________________________ Координати стаціонарної точки № 3: x1=0 x2=0 x3=0 Похідна поля напрямків: df1 / dx1 = -10 df1 / dx2 = 10 df1 / dx3 = 0 df2 / dx1 = r df2 / dx2 = -1 df2 / dx3 = 0 df3 / dx1 = 0 df3 / dx2 = 0 df3 / dx3 = ( - 8)/3 Характеристичне рівняння: (3*z**3 + 41*z**2 + 2*z*( - 15*r + 59) + 80*( - r + 1))/3 = 0 Власні значення - корені характеристичного рівняння: z=(sqrt(40*r + 81) - 11)/2 z=( - sqrt(40*r + 81) - 11)/2 z=( - 8)/3
При r < 1 система Лоренца (1) має єдину стаціонарну (нерухому) точку — початок координат, до якого притягуються всі траекторії (рис. 4).
При r > 1 початок координат вже нестійкий, бо відповідна матриця ∂ f (x) / ∂ x має одне додатне і два від’ємних власних числа. Крім початку координат при r > 1 існує ще дві стаціонарні точки:
Спочатку всі власні числа λ матриці ∂ f / ∂ x у двох нових стаціонарних точках від’ємні (рис. 5).
Згодом при зростанні r пара від’ємних власних чисел перетворюється на пару комплексно спряжених чисел з від’ємною дійсною частиною. Кожна з двох траекторії G1 і G2, що «виходять» з початку координат у напрямках власного вектора матриці ∂ f (0) / ∂ x з додатним власним числом, і «входили» у стійку стаціонарну точку X1 чи X2, починають закручуватися спіраллю навколо цієї стаціонарної точки (рис. 6).
При r ≈ 13,92 траекторії Г1 і Г2 «входять» у початок координат у напрямку у напрямку власного вектора матриці ∂ f (0) / ∂ x з від’ємним власним числом (рис. 7).
При подальшому зростанні r в околі ще стійких стаціонарних точок з’являються два нестійкі цикли Ф1 і Ф2. Лінійні частини операторів слідування для цих циклів мають по два власних значення, що лежать по різні боки від 1. Інакше кажучи, по одному напряму тракторії притягуються до цих циклів, а по іншому — відштовхуються. При цьому дві траекторії G1 і G2 притягуються вже до інших стаціонарних точок, ніж до появи нестійких циклів, і закручуються спіраллю навколо них (рис. 8).
При r ≈ 24,06 траекторії G1 і G2 потрапляють на многовиди притягування до нестійких циклів Ф1 і Ф2.
При r = r0 = σ(σ + b + 3)/(σ – b – 1) ≈ 24,74 матриця ∂ f (Xj) / ∂ x (j = 1, 2) має два комплексно спряжених власних числа на уявній осі. Дійсні частини цих власних чисел зростають при зростанні r. Стаціонарні точки X1 та X2 гублять свою стійкість, поглинаючи нестійкі цикли (так звана біфуркація Пуанкаре — Андронова — Хопфа).
Вже при r > 13.92 система Лоренца має граничну інваріантну множину, яка при r < r0 не є стійкою, тобто не притягує траекторії. При r ∈ (r0; 50]
ця множина стає стійкою. Її називають аттрактором Лоренца. На рис. 9 подано двовимірну проекцію траекторії в околі аттрактора Лоренца. Траекторія робить кілька обертів навколо однієї стаціонарної точки, потім — навколо іншої і т.і. Кількість таких послідовних обертів навколо стаціонарної точки настільки істотно залежить від координат початкової точки, що виникає враження хаотичності руху і, навіть, випадковості такого руху.
Аттрактор Лоренца A має такі властивості (обґрунтування цих властивостей полягає у тлумаченні результатів числових розрахунків).
A є аттрактором у тому розумінні, що існує відкрита в R3 множина B, при якій
| tA = | ∩ | TtB |
| t ≥ 0 |
Тут Tt — оператор зсуву вздовж траекторій системи на час t. Інакше кажучи, всі траекторії, що мають початок у B (для системи Лоренца за B можна взяти все R3 без стаціонарних точок), притягуються до A.
A містить всюди щільну множину періодичних траекторій, кожна з яких є нестійкою.
Для довільної траекторії в A з трьох різних характиристичних показників Ляпунова один дорівнює 0 (для напряму потоку f ), один — від’ємний (вздовж відповідного напряму відбувається притягування траекторій), а ще один — додатний (вздовж відповідного напряму відбувається експотенціалльне «розбігання» траекторій навіть при як завгодно малих відхиленнях початкових значень).
Локально аттрактор Лоренца має структуру добутку інтервалу на досконалу канторову множину.
При більших величинах r виявлено стійкі періодичні розв'язки. Наприклад, при r = 234 система Лоренца має стійкий цикл. При зменшенні величини r такий цикл втрачає стійкість, породжуючи стійкий цикл подвоєного періоду. При подальшому зменшенні величини r новостворений цикл втрачає стійкість, породжуючи новий стійкий цикл подвоєного періоду і т.д. і т.і. Таким чином виникає нескінчена спадна послідовність {rk} величин параметра, при якому система Лоренца демонструє (біфуркацію) подвоєння періоду стійкого циклу. Виявляється, ця послідовність задовольняє закону універсальності Фейгенбаума, тобто рівності:
|
Тут δ = 4.6692… — певна універсальна стала, для якої доведено відповідне твердження для неперервних перетворень прямої (інтервала) в себе (Вул Е. Б., Синай Я. Г., Ханин К. М., Универсальность Фейгенбаума и термодинамический формализм // УМН, 1984, 39:3(237), с. 3–37).