В первой части мы обсудили программу для моделирования логистического отображения и основные режимы поведения, которые может демонстрировать эта проста модель динамической системы. Теперь мы пойдем дальше и немного усовершенствуем программу, созданную до этого. Дело в том, что там мы просто численно моделировали численность популяции с заранее заданными временем установления и временем рассмотрения. Такой подход создаст нам трудности в будущем, поскольку время установления может быть очень разным. Чем ближе к точке бифуркации, тем время переходного процесса больше. Обсудим в этой статье модификацию программы и зависимости времени переходного процесса от управляющего параметра при приближении к точке бифуркации.
Во-первых немного поменяем смысл параметра transient_time, пусть это будет не время переходного процесса, а максимально допустимое время переходного процесса max_transient. Также добавим автоматическое увеличение параметра observation_time в случаях очень большого периода колебания и хаотических колебаний. Программа делает вывод только о периодичности колебаний, если ей хватает времени transient_time, если времени не хватило, то вывод не делается, потому что это либо хаотические колебания, либо точка бифуркации, либо просто не хватило времени обнаружить периодичность в колебаниях. Код программы приведен в конце статьи. В комментарии в начале программы приведены примеры компиляции и запуска программы. Компилируется программа командой "g++ -o AdapLog AdapLog.cpp -lm", а запускается командой "./AdapLog 2.999999 0.5 100000000 100", где первое число задает параметр отображения, второе определяет начальное состояние, третье ограничивает время переходного процесса сверху, а последнее задает количество точек в конце процесса, которое выведется в файл в конце вычислений. В дополнение к временной реализации AdapLog также выводит информацию о полученном режиме, если он детектирован как периодический.
Выполним в консоли bash команду "for r in $(LC_ALL=C seq 0.0001 0.0001 4.00); do ./AdapLog ${r} 0.5 1000000 10; cat time-series.dat >> time-series-all.dat; done". Последняя команда будет очень долго выполняться, несколько дней, но можно выбирать либо меньшие отрезки (например, "seq 3.700001 0.000001 3.80"), либо более грубый шаг (например, "seq 0.01 0.01 4.00"). Скрипт сгенерирует очень большой файл, который лучше немного почистить перед построением в gnuplot (скрипты для построения всех графиков, представленных в этой статье, приведены в конце статьи).
Зависимость времени установления от близости к точке бифуркации
На рисунке 1 приведена зависимость времени установления режима от близости к точке бифуркации, которая иллюстрирует тот факт, что, при приближении к точкам бифуркации, меняется характер устойчивости режимов. Каждый раз перед удвоением периода происходит бифуркация, ровно в которой (если параметр r задать равным точно бифуркационному значению) существуют оба режима (тот, что реализовывался перед бифуркацией, и тот, который будет наблюдаться после бифуркации). При этом оба режима характеризуются нейтральной устойчивостью, поэтому время выхода на любой из режимов стремится к бесконечности. Вблизи же точек бифуркации время выхода на наблюдаемый режим значительно возрастает, но является конечным. Интересно отметить, что в точке r=1 не происходит удвоения периода, а возникает цикл периода 1. Это тоже бифуркация, поэтому время выхода на режим вблизи r=1 также резко увеличивается.
Бифуркационная диаграмма
Бифуркационная диаграмма на рисунке 2 представляет собой набор точек временной реализации установившегося режима для различных значений r, отложенных по горизонтали. На диаграмме хорошо видно, что от r=0 до 1 устанавливается режим со значениями x(n)=0. При r>1 появляется цикл периода 1 и увеличивается равновесное значение переменной x. При r=3 происходит удвоение периода, которое затем повторяется до появления режима динамического хаоса. Стоит также заметить, что после достижения режима динамического хаоса вновь могут наблюдаться периодические циклы разных периодов. Строго говоря, в этой области хаоса также присутствуют все циклы любых периодов из порядка Шарковского (ну, то есть отвечающие всем натуральным числам).
Зависимость периода от управляющего параметра
Представленная на рисунке 3 зависимость по смыслу сильно схожа с бифуркационной диаграммой, только по вертикали либо отложен период колебаний, либо точка отсутствует, если период не обнаружен. Само изображение анализировать сложно, но исходный файл данных можно просмотреть и найти вручную все обнаруженные бифуркационные переходы.
В первом каскаде бифуркаций удвоения периода начинаются с цикла периода 1 и происходят при r=3, 3.44949, 3.54409, 3.564407, 3.568759, 3.5696915, 3.569892 и т.д. В этом ряду отношения Фейгенбаума равны 4.7515, 4.6562, 4.6684, 4.6670 и 4.6509. По первым трем элементам получившегося ряда видно, что отношение стремится к константе Фейгенбаума (4.6692), после чего (начиная с 4-го элемента) ошибка возрастает, что связано с ошибкой численного алгоритма, который использовался для определения периода колебаний.
Особенностью хаотических систем со сценарием Фейгенбаума перехода к хаосу является то, что, после выхода на хаотический режим, при продолжении увеличения параметра, можно вновь обнаружить периодические колебания самых разных периодов. Такие диапазоны, где наблюдаются периодические колебания, называют окнами периодичности. Окна периодичности могут быть как очень узкими, так и достаточно широкими. Самое широкое окно - окно периодичености 3. В каждом таком окне вновь наблюдается переход к хаосу через каскад бифуркаций удвоения периода. В рассматриваемой нами модели наблюдаются следующие окна периодичности (обнаружены при проходе с шагом 0.0001 по параметру r):
На основе периода 24: 3.572533⪅r⪅3.572800.
На основе периода 20: 3.577513⪅r⪅3.577753.
На основе периода 12: 3.582024⪅r⪅3.583285.
На основе периода 10: 3.605209⪅r⪅3.606366 и 3.647007⪅r⪅3.647279.
На основе периода 6: 3.626554⪅r⪅3.632714.
На основе периода 8: 3.662109⪅r⪅3.662652 и 3.800741⪅r⪅3.800943.
На основе периода 9: 3.687197⪅r⪅3.687322; 3.717094⪅r⪅3.717244 и 3.853614⪅r⪅3.853986.
На основе периода 7: 3.701641⪅r⪅3.702485 и 3.774134⪅r⪅3.774658.
На основе периода 5: 3.738173⪅r⪅3.743001 и 3.905572⪅r⪅3.906449.
На основе периода 3: 3.828428⪅r⪅3.849429.
На основе периода 4: 3.960102⪅r⪅3.961193.
Если бы в поиске окон периодичности мы бы прошли с меньшим шагом, то обнаружили бы еще больше окон периодичности, потому что характер изменения динамики хаотической системы представляется фрактальной структурой.
Заключение
Представленный способ исследования сценария перехода к хаосу через каскад бифуркаций удвоения периода не является самым изящным, но он довольно простой. Для обнаружения бифуркаций (более достоверно) нужно считать мультипликаторы или показатели Ляпунова, а не просто период колебаний.
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
Программа AdapLog для моделирования поведения логистического отображения с адаптивным временем установления
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <cmath>
#include <iomanip>
/*
/*The program written by Andrei Bukh, to compile it run:
* $ g++ -o AdapLog AdapLog.cpp -lm
* $ ./AdapLog 2.999999 0.5 100000000 100
*/
int logistic (double *in, double *out, double par_r) {
if (((*in) < 0.0) || ((*in) > 1.0)) return 1;
if ((par_r < 0.0) || (par_r > 4.0)) return 2;
(*out) = par_r * (*in) * (1.0 - (*in));
return 0;
}
int main (int argc, char* argv[]) {
double parameter_r = 0.0;
if (argc > 1) {
bool rightnumber = true;
try {
parameter_r = std::stod (argv[1]);
} catch (std::invalid_argument) {
rightnumber = false;
}
if (!rightnumber) {
std::cout << "Please, enter parameter r = ";
std::cin >> parameter_r;
}
}
else {
std::cout << "Please, enter parameter r = ";
std::cin >> parameter_r;
}
int maxdelay = 1000;
double *x; x = new double [maxdelay];
x[0] = 0.0;
if (argc > 2) {
bool rightnumber = true;
try {
x[0] = std::stod (argv[2]);
} catch (std::invalid_argument) {
rightnumber = false;
}
if (!rightnumber) {
std::cout << "Please, enter x_0 = ";
std::cin >> x[0];
}
}
else {
std::cout << "Please, enter x_0 = ";
std::cin >> x[0];
}
long int max_transient = 0;
if (argc > 3) {
bool rightnumber = true;
try {
max_transient = std::stoi (argv[3]);
} catch (std::invalid_argument) {
rightnumber = false;
}
if (!rightnumber) {
std::cout << "Please, enter max_transient = ";
std::cin >> max_transient;
}
}
else {
std::cout << "Please, enter max_transient = ";
std::cin >> max_transient;
}
int observation_time = 0;
if (argc > 4) {
bool rightnumber = true;
try {
observation_time = std::stoi (argv[4]);
} catch (std::invalid_argument) {
rightnumber = false;
}
if (!rightnumber) {
std::cout << "Please, enter observation_time = ";
std::cin >> observation_time;
}
}
else {
std::cout << "Please, enter observation_time = ";
std::cin >> observation_time;
}
std::cout << "parameter_r = " << std::setprecision(12) << parameter_r << "\n";
std::cout << "x_0 = " << std::setprecision(12) << x[0] << "\n";
std::cout << "max_transient = " << max_transient << "\n";
std::cout << "observation_time = " << observation_time << "\n";
int marker = 0;
long int time = 0;
long int i = 0;
bool periodicityisnotfound = true;
double glob_min = 1.0e6, glob_max = -1.0e6;
int glob_period = 0, suc_counter = 0;
while ((i<max_transient)&&(periodicityisnotfound)) {
double min = 1.0e6;
double max = -1.0e6;
int startnum = marker+1;
if (startnum==maxdelay) startnum=0;
for (int j=0;j<maxdelay;j++) {
int pos = startnum + j; if (pos>=maxdelay) pos -= maxdelay;
if (x[pos]<min) min = x[pos];
if (x[pos]>max) max = x[pos];
}
double nearestx = 1.0e6;
int period = 0;
for (int j=1;j<maxdelay;j++) {
int pos = startnum + j; if (pos>=maxdelay) pos -= maxdelay;
if ((fabs(x[startnum]-x[pos])+1.0e-6)<(fabs(x[startnum]-nearestx))) {
nearestx = x[pos];
period = j;
}
}
if (
(fabs(min-glob_min)>1.0e-20) ||
(fabs(max-glob_max)>1.0e-20) ||
(period!=glob_period)
) {
glob_min = min;
glob_max = max;
glob_period = period;
suc_counter = 0;
}
else {
suc_counter++;
static int time_moment = 0;
if (suc_counter==1) time_moment = time - maxdelay;
if (suc_counter > (1000)) {
std::ofstream PFF("transient-time.dat", std::ios_base::app);
PFF << std::setprecision(12) << parameter_r << " " << glob_period << " " << time_moment << " " << glob_min << " " << glob_max << "\n";
PFF.close();
periodicityisnotfound = false;
std::cout << "The cycle of period " << glob_period << " is observed\n";
}
}
if ((time%1000000)==0)
std::cout << "time=" << time << "\n";
int nextmarker = marker + 1; if (nextmarker==maxdelay) nextmarker = 0;
int error = logistic (&x[marker], &x[nextmarker], parameter_r);
if (error == 1) {
std::cout << "Variable x has gone to infinity\n";
return error;
}
if (error == 2) {
std::cout << "Parameter r has to be in range [0:4]\n";
return error;
}
i++;
time++;
marker = nextmarker;
}
std::ofstream TS_file("time-series.dat");
if (glob_period>observation_time) observation_time = glob_period;
if (periodicityisnotfound) if (observation_time<100) observation_time = 100;
for (int i=0; i<observation_time; i++) {
int nextmarker = marker + 1; if (nextmarker==maxdelay) nextmarker = 0;
int error = logistic (&x[marker], &x[nextmarker], parameter_r);
if (error == 1) {
std::cout << "Variable x has gone to infinity\n";
return error;
}
if (error == 2) {
std::cout << "Parameter r has to be in range [0:4]\n";
return error;
}
TS_file << parameter_r << " " << time << " " << std::setprecision(12) << x[marker] << "\n";
time++;
marker = nextmarker;
}
TS_file << parameter_r << " " << time << " " << std::setprecision(12) << x[marker] << "\n";
TS_file.close();
delete [] x;
return 0;
}
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
Скрипты для построения графиков этой статьи в gnuplot
Для сортировки и очистки файлов в bash:
sort file1.dat | uniq -u > file2.dat
LC_ALL=C sort -g -k3,3 -k1,1 file1.dat > file2.dat
Для построения изображений:
set autoscale; unset label; unset log y; set xrange [3.82:3.87]; unset key; set lmargin screen 0.01; set bmargin screen 0.02; set rmargin screen 0.99; set tmargin screen 0.96; unset xlabel; unset ylabel; unset tics; unset border; plot "time-series-all.dat" u 1:3 w p pt 5 ps 0.2 lc rgb "forest-green" title "Variable values vs parameter r for logistic map"
set autoscale; unset label; set border 15; set log y; set label "Transient time vs parameter r for logistic map" at screen 0.143,0.85 left font "Helvetica, 30"; set lmargin screen 0.08; set bmargin screen 0.2; set rmargin screen 0.97; set tmargin screen 0.94; set xlabel "r" offset 0,-1.5 font "Helvetica, 30"; set ylabel "t_{tr}" offset -6.5,0 font "Helvetica, 30"; set tics font "Helvetica, 30"; set xtics 0.5 offset 0,-0.5; set ytics 100; plot "transient-time.dat" u 1:3 w p pt 7 ps 0.4 lc rgb "forest-green" notitle
set autoscale; unset label; set border 15; set log y; set label "Transient time vs parameter r for logistic map" at screen 0.143,0.85 left font "Helvetica, 30"; set lmargin screen 0.08; set bmargin screen 0.2; set rmargin screen 0.97; set tmargin screen 0.94; set xlabel "r" offset 0,-1.5 font "Helvetica, 30"; set ylabel "t_{tr}" offset -6.5,0 font "Helvetica, 30"; set tics font "Helvetica, 30"; set xtics 0.5 offset 0,-0.5; set ytics 100; plot "transient-time2.dat" u 1:3 w p pt 7 ps 0.6 lc rgb "forest-green" notitle
set autoscale; set xrange [2.95:4]; unset label; set border 15; set log y; set label "Period vs parameter r for logistic map" at screen 0.143,0.85 left font "Helvetica, 30"; set lmargin screen 0.08; set bmargin screen 0.2; set rmargin screen 0.97; set tmargin screen 0.94; set xlabel "r" offset 0,-1.5 font "Helvetica, 30"; set ylabel "T" offset -6.5,0 font "Helvetica, 30"; set tics font "Helvetica, 30"; set xtics 0.1 offset 0,-0.5; set ytics 4; plot "transient-time2.dat" u 1:2 w p pt 7 ps 0.6 lc rgb "forest-green" notitle
set autoscale; set xrange [2.95:3.7]; unset label; set border 15; unset log y; set label "Variable values vs parameter r for logistic map" at screen 0.143,0.85 left font "Helvetica, 30"; set lmargin screen 0.08; set bmargin screen 0.2; set rmargin screen 0.97; set tmargin screen 0.94; set xlabel "r" offset 0,-1.5 font "Helvetica, 30"; set ylabel "x" offset -6.5,0 font "Helvetica, 30"; set tics font "Helvetica, 30"; set xtics 0.2 offset 0,-0.5; set ytics 0.2; plot "time-series-all2.dat" u 1:3 w p pt 7 ps 0.4 lc rgb "forest-green" notitle
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.