В нелинейной динамике есть несколько способов определения устойчивости наблюдаемых режимов. Фазовая траектория устойчива по Лагранжу, если можно выделить конечную область фазового пространства, которую рассматриваемая траектория не покидает. В этом смысле хаотические аттракторы могут быть устойчивыми, если, несмотря на нерегулярность поведения, они целиком заключены в некоторый фазовый объем конечных размеров. Устойчивость по Пуассону тесно связана с теорией возвратов Пуанкаре. Если фазовая траектория время от времени возвращается в ближайшие окрестности своих предыдущих состояний, она устойчива по Пуассону. Сложнее этих двух методов определяется устойчивость по Ляпунову. Траектория в фазовом пространстве является устойчивой в смысле Ляпунова, если небольшие возмущения в фазовые состояния вдоль траектории уменьшаются с течением времени. Другими словами, траектория устойчива по Ляпунову, если возмущенные траекотрии со временем сливаются с невозмущенной траекторией.
На практике устойчивость траекторий можно определить двумя методами. Первый вариант подразумевает рассмотрение оригинальной динамической системы и ее копий. При этом в фазовые состояния копий динамической системы вводят отклонения, а затем наблюдают, как ведут себя копии по отношению к оригиналу. Если отклонения уменьшаются, фазовая траектория в оригинальной системе является устойчивой. Этот метод требует настройки для каждой динамической системы, поскольку непонятно, например, что считать малым отклонением. Второй способ более изящен, но для этого нужно рассмотреть оригинальную динамическую систему с линеаризованными копиями. Тогда, если в линеаризованном фазовом пространстве ненулевой вектор уменьшается, фазовая траектория является устойчивой по Ляпунову. В этом случае важно только то, чтобы начальный вектор был ненулевым, его величина в линеаризованной системе не имеет значения. Второй случай, кроме того, удобен тем, что в нем относительно просто найти спектр показателей Ляпунова. Тем, кому интересен вопрос, рекомендую книгу С.П. Кузнецова "Динамический хаос".
Рассмотрим устойчивость на примере генератора хаоса Анищенко-Астахова. Это трехмерная автономная динамическая система, так что для нее можно определить максимум три показателя Ляпунова. Программирую я на C/C++, поэтому в начале подключим базовые библиотеки C++.
#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
Определим количество независимых переменных (размер фазового пространства) и общее количество переменных (основное фазовое пространство + 3 линеаризованных пространства для вычисления показателей Ляпунова).
int systemvars = 3;
int totalvars = systemvars*(systemvars+1);
int totalpars = 2;
Уравнения генератора Анищенко-Астахова вместе с соответствующими уравнениями в вариациях
представлены функцией int func (...).
int func (double *** inout_vars, int mkr, int srfr, int stto, double * in_times, double * in_pars) {
double x = inout_vars[mkr][srfr][0];
double ax = fabs(x);
double y = inout_vars[mkr][srfr][1];
double z = inout_vars[mkr][srfr][2];
double m = in_pars[0];
double g = in_pars[1];
inout_vars[mkr][stto][0] = x*(m-z) + y;
inout_vars[mkr][stto][1] = - x;
inout_vars[mkr][stto][2] = g*(0.5*x*(x+ax)-z);
for (int v=systemvars;v<totalvars;v+=systemvars) {
double tx = inout_vars[mkr][srfr][v+0];
double ty = inout_vars[mkr][srfr][v+1];
double tz = inout_vars[mkr][srfr][v+2];
inout_vars[mkr][stto][v+0] = tx*(m-z) + ty - tz*x;
inout_vars[mkr][stto][v+1] = - tx;
inout_vars[mkr][stto][v+2] = g*(tx*(x + ax) - tz);
}
return 0;
}
Функция int func (...) используется методом интегрирования (4-шаговый метод Рунге—Кутты) int rk4step (...).
int rk4step (double *** inout_vars, double *times, double *pars, int marker, int md, double h) {
int pm = marker-1; if (pm<0) pm+=md;
for (int v=0;v<totalvars;v++)
inout_vars[marker][5][v] = inout_vars[marker][0][v] = inout_vars[pm][5][v];
times[marker] = times[pm]; func (inout_vars,marker,0,4,times,pars); /*K1*/
for (int v=0;v<totalvars;v++)
inout_vars[marker][5][v] += h * inout_vars[marker][4][v] / 6.0;
for (int v=0;v<totalvars;v++)
inout_vars[marker][1][v] = inout_vars[marker][0][v] + h * 0.5 * inout_vars[marker][4][v];
times[marker] += 0.5*h; func (inout_vars,marker,1,4,times,pars); /*K2*/
for (int v=0;v<totalvars;v++)
inout_vars[marker][5][v] += h * inout_vars[marker][4][v] / 3.0;
for (int v=0;v<totalvars;v++)
inout_vars[marker][2][v] = inout_vars[marker][0][v] + h * 0.5 * inout_vars[marker][4][v];
func (inout_vars,marker,2,4,times,pars); /*K3*/
for (int v=0;v<totalvars;v++)
inout_vars[marker][5][v] += h * inout_vars[marker][4][v] / 3.0;
for (int v=0;v<totalvars;v++)
inout_vars[marker][3][v] = inout_vars[marker][0][v] + h * inout_vars[marker][4][v];
times[marker] += 0.5*h; func (inout_vars,marker,3,4,times,pars); /*K4*/
for (int v=0;v<totalvars;v++)
inout_vars[marker][5][v] += h * inout_vars[marker][4][v] / 6.0;
return 0;
}
Метод Рунге—Кутты нужен для численного моделирования генератора Анищенко-Астахова. Если в эту функцию передать массив с переменными и номерами входного и выходного подмассивов, будет выполнен 1 шаг интегрирования методом Рунге—Кутты. Функция int init_state (...) предназначена для инициализации фазового состояния и начальных векторов возмущения
для вычисления показателей Ляпунова. В этом примере начальные вектора заданы в виде единичных базовых векторов фазового пространства.
Можно задать их и иначе, но важно, чтобы все они были попарно ортогональны.
int init_state (double *times, double ***outvars, double *invars, double *outlyaps, double *outpars, double *inpars, int md) {
times[0] = 0.0;
for (int d=0;d<md;d++)
for (int s=0;s<6;s++)
for (int v=0;v<totalvars;v++)
if (v<systemvars) {
outvars[d][s][v] = invars[v];
outlyaps[v] = 0.0;
}
else
if ((v+1)%(systemvars+1)==0)
outvars[d][s][v] = 1.0;
else
outvars[d][s][v] = 0.0;
for (int p=0;p<totalpars;p++) outpars[p] = inpars[p];
return 0;
}
Вспомогательная функция int calc_lengths (...) нужна для вычисления длин векторов возмущения.
int calc_lengths (double *lng, double ***in_vars, int mrk) {
for (int v=0;v<systemvars;v++) {
lng[v]=0.0;
for (int v_position=0;v_position<systemvars;v_position++) {
lng[v]+=in_vars[mrk][0][systemvars*(1+v)+v_position]*in_vars[mrk][0][systemvars*(1+v)+v_position];
}
lng[v]=sqrt(lng[v]);
}
return 0;
}
Функция int calc_angles (...) вычисляет углы между векторами возмущения, чтобы контролировать их неколлинеарность.
Дело в том, что коллинеарные векторы нельзя оргонализовать, а для вычисления младших показателей Ляпунова нужна ортогонализация векторов.
int calc_angles (double *out_angles, double *out_prods, double *in_lngths, double ***vars, int mrk) {
for (int v=0;v<systemvars;v++) {
out_prods[v]=0.0;
for (int v_position=0;v_position<systemvars;v_position++) {
if (v<(systemvars-1)) out_prods[v]+=vars[mrk][0][systemvars*(1+v)+v_position]*vars[mrk][0][systemvars*(2+v)+v_position];
}
}
for (int v=0;v<systemvars;v++) {
out_angles[v]=0.0;
if (v<(systemvars-1)) out_angles[v]=acos(out_prods[v]/(in_lngths[v]*in_lngths[v+1]));
}
return 0;
}
Ключевой функцией для вычисления спектра показателей Ляпунова является int gramshmidt (...), которая выполняет ортогонализацию и перенормировку векторов.
Без этой процедуры (или аналога) вектора возмущения выстоятся в направлении максимального роста возмущения и дадут три одинаковых значения.
int gramshidt (double ***outvars, double *out_prods, double *out_lng, double *out_lyaps, int md, int mrk) {
for (int v=0;v<systemvars;v++) {
for (int d=0;d<md;d++)
for (int s=0;s<6;s++) {
for (int v_subtructed=0;v_subtructed<v;v_subtructed++) {
out_prods[v_subtructed]=0.0;
for (int v_position=0;v_position<systemvars;v_position++) {
out_prods[v_subtructed]+=outvars[d][s][systemvars*(1+v)+v_position]*outvars[d][s][systemvars*(1+v_subtructed)+v_position];
}
}
for (int v_subtructed=0;v_subtructed<v;v_subtructed++) {
for (int v_position=0;v_position<systemvars;v_position++) {
outvars[d][s][systemvars*(1+v)+v_position]-=out_prods[v_subtructed]*outvars[d][s][systemvars*(1+v_subtructed)+v_position];
}
}
}
calc_lengths(out_lng, outvars, mrk);
out_lyaps[v] += log(out_lng[v]);
for (int d=0;d<md;d++)
for (int s=0;s<6;s++)
for (int v_position=0;v_position<systemvars;v_position++)
outvars[d][s][systemvars*(1+v)+v_position]/=out_lng[v];
out_lng[v]=0.0;
for (int v_position=0;v_position<systemvars;v_position++) {
out_lng[v]+=outvars[mrk][0][systemvars*(1+v)+v_position]*outvars[mrk][0][systemvars*(1+v)+v_position];
}
out_lng[v]=sqrt(out_lng[v]);
}
return 0;
}
В функции int singleparset (...) выполняются все основные вычисления, подготавливаются массивы данных, вызывается метод интегрирования, контролируются ортогонализации и перенормировки векторов возмущения. С целью обеспечения хорошей точности вектора возмущения интегрируются с самого начала (включая переходный процесс), а сразу после времени установления вектора ортогонализуются и нормируются в очередной раз, с которого начинает собираться сумма логарифмов в показатели Ляпунова.
int singleparset (double *trvar, double *trpar, int maxdelay, double time, double trans, double timestep) {
double*** vars; double* times; double* pars; double *lyaps, *lngths, *prods, *angles;
vars = new double** [maxdelay];
for (int d=0;d<maxdelay;d++) vars[d] = new double* [6];
for (int d=0;d<maxdelay;d++)
for (int s=0;s<6;s++) vars[d][s] = new double [totalvars];
times = new double [maxdelay];
pars = new double [totalpars];
lyaps = new double [systemvars];
lngths = new double [systemvars];
prods = new double [systemvars];
angles = new double [systemvars];
init_state (times, vars, trvar, lyaps, pars, trpar, maxdelay);
std::ofstream ts_file("TS.dat"); ts_file << std::setprecision(9) << std::scientific;
int marker = 0;
long t = 0;
double ct = .0;
bool lyap_started = false; double lyap_start_time = 0.0; bool calconemoretime = true; bool calctwotimesmore = true;
while (calconemoretime) {
ct += timestep; if (!(ct < time)) { if (calctwotimesmore) calctwotimesmore = false; else calconemoretime = false; }
t++;
marker++; if (marker>=maxdelay) marker = 0;
int prevmarker = marker-1; if (prevmarker<0) prevmarker = maxdelay-1;
if (((t%(1000000))==0)) std::cout<<"time="<<t*timestep<<"\n";
rk4step (vars, times, pars, marker, maxdelay, timestep);
if (ct>trans) {
if (!lyap_started) {
lyap_started = true;
lyap_start_time = times[marker];
gramshidt(vars, prods, lngths, lyaps, maxdelay, marker);
for (int v=0;v<systemvars;v++) lyaps[v] = 0.0;
}
else {
bool rebuild = false;
if (!calctwotimesmore) rebuild = true;
calc_lengths(lngths, vars, marker);
for (int v=0;v<systemvars;v++)
if ((fabs(lngths[v])>10.0)||(fabs(lngths[v])<0.1)) rebuild = true;
calc_angles(angles, prods, lngths, vars, marker);
for (int v=0;v<systemvars;v++)
if ((angles[v]<1.39626)||(angles[v]>1.74533)) rebuild = true;
if (rebuild)
gramshidt(vars, prods, lngths, lyaps, maxdelay, marker);
}
ts_file << pars[0] << " " << pars[1] << " " << times[marker];
for (int i=0; i<systemvars; i++) ts_file << " " << vars[marker][0][i];
for (int i=0; i<systemvars; i++) ts_file << " " << (lyaps[i]+log(lngths[i]))/(times[marker]-lyap_start_time);
ts_file << "\n";
}
else {
gramshidt(vars, prods, lngths, lyaps, maxdelay, marker);
for (int v=0;v<systemvars;v++) lyaps[v] = 0.0;
}
}
ts_file.close();
std::ofstream lyapf ("lyaps.dat");
lyapf << pars[0] << " " << pars[1] << " " << times[marker];
for (int i=0; i<systemvars; i++) lyapf << " " << vars[marker][0][i];
for (int i=0; i<systemvars; i++) lyapf << " " << (lyaps[i]+log(lngths[i]))/(times[marker]-lyap_start_time);
lyapf << "\n";
lyapf.close();
delete [] angles;
delete [] prods;
delete [] lngths;
delete [] lyaps;
delete [] pars;
delete [] times;
for (int i=0;i<maxdelay;i++)
for (int j=0;j<6;j++) delete [] vars[i][j];
for (int i=0;i<maxdelay;i++) delete [] vars[i];
delete [] vars;
return 0;
}
Две следующие функции выполняют вспомогательную функцию программы и считывают аргументы командной строки для передачи их в динамическую систему в качестве параметров.
float takefloatpar (int argc, char* argv[], int num) {
double res = 0.0;
if (argc > num) {
bool rightnumber = true;
try {
res = std::stod (argv[num]);
} catch (std::invalid_argument) {
rightnumber = false;
}
if (!rightnumber) {
std::cout << "Please, enter parameter " << num << ": ";
std::cin >> res;
}
}
else {
std::cout << "Please, enter parameter " << num << ": ";
std::cin >> res;
}
std::cout << num << " parameter is " << res << "\n";
return res;
}
long int takeintpar (int argc, char* argv[], int num) {
long int res = 0;
if (argc > num) {
bool rightnumber = true;
try {
res = std::stoi (argv[num]);
} catch (std::invalid_argument) {
rightnumber = false;
}
if (!rightnumber) {
std::cout << "Please, enter parameter " << num << ": ";
std::cin >> res;
}
}
else {
std::cout << "Please, enter parameter " << num << ": ";
std::cin >> res;
}
std::cout << num << " parameter is " << res << "\n";
return res;
}
Наконец, функция int main (...) принимает аргументы командной строки, парсит их и передает в основную функцию программы.
int main (int argc, char* argv[]) {
double pars[totalpars];
pars[0] = takefloatpar(argc, argv, 1);
pars[1] = takefloatpar(argc, argv, 2);
double vars[totalvars];
vars[0] = takefloatpar(argc, argv, 3);
vars[1] = takefloatpar(argc, argv, 4);
vars[2] = takefloatpar(argc, argv, 5);
double time = takefloatpar(argc, argv, 6);
double trans = takefloatpar(argc, argv, 7);
double timestep = takefloatpar(argc, argv, 8);
singleparset(vars,pars,2,time,trans,timestep);
return 0;
}
Запустить программу можно, например, выполнив команду ./GIN-LCE 1.563 0.17 0.1 0.1 0.1 2000 1000 0.001. Тогда будет выполнено моделирование генератора Анищенко-Астахова с параметрами m=1.563 и g=0.17. Начальными условиями будут x0=0.1, y0=0.1, z0=0.1. Программа будет производить вычисления в течение 2000 единиц безразмерного времени и отбросит при этом 1000 единиц. Шаг интегрирования установлен равным 0.001.
На рисунке 1 приведен результат вычисления показателей Ляпунова для хаотического аттрактора в генераторе Анищенко-Астахова при m=1.563 и g=0.17. Для диагностики режима поведения генератора нужны установившиеся значения показателей (значения при t=2000 в данном случае), которые равны значениям 0.03196, 0.0000395 и -0.166. Первый показатель Ляпунова принимает положительное значение и подтверждает хаотичность установленного режима. Второй показатель Ляпунова можно считать нулевым (когда реализуется решение, отличное от состояния равновесия, хотя бы один показатель Ляпунова равен нулю, потому что направление вдоль траектории является нейтрально устойчивым). Третий показатель Ляпунова явялется отрицательным и свидетельствует о том, что в фазовом пространстве есть диссипация, которая приводит к возвращению траектории в область хаотического аттрактора.
Для исследования поведения модели с помощью показателей Ляпунова достаточно запустить скрипт, который будет вызывать программу и сохранять значения показателей Ляпунова в отдельный файл.
#!/bin/bash
for m in $(LC_ALL=C seq 1.4000 0.0050 1.4000); do
for g in $(LC_ALL=C seq 0.1000 0.0005 0.2000); do
./GIN $m $g 0.1 0.1 0.1 2000 1000 0.001;
cat lyaps.dat >> l_all.dat
done;
done;
exit
p.s. в представленной здесь программе, в сравнении с алгоритмом, взятым с книги С.П. Кузнецова "Динамический хаос", добавлен контроль угла между векторами, так что, при отклонении на 10 градусов от прямого угла между векторами возмущения, происходит внеочередная ортогонализация и перенормировка. Основная перенормировка происходит при значительном удлинении или сокращении одного из векторов. Такой алгоритм, на мой взгляд, гораздо более гибок.
Подписывайтесь на мой канал, эта статья не из простых, но здесь планирую добавлять очень разный материал по физике и математике. Задавайте вопросы в комментариях и пишите о своих пожеланиях.