Здесь были рассмотрены изменения в спектре, происходящие при реализации перехода к хаосу через каскад бифуркаций удвоения периода на примере генератора с инерционной нелинейностью Анищенко-Астахова. Но это совсем не единственная модель, которая такой сценарий иллюстрирует. Многообразие моделей, показывающих один и тот же феномен, является всегда косвенным признаком важности наблюдаемого явления. При этом поведение различных моделей количественно различается, но отдельные качественные признаки повторяются. Рассмотрим в этой статье генератор Ресслера. Это модель, которая появилась совсем не в результате попыток моделирования реальных динамических систем, как например, генератор с инерционной нелинейностью, рассмотренный в предыдущей статье. Однако модель Ресслера проще и поэтому заслуживает особого внимания. Одной из задач нелинейной динамики, как междисциплинарной науки, является выявление принципиальных особенностей моделей, реализующих один и тот же феномен. Поэтому, чем проще модель, в которой возможно наблюдение некоторого феномена, тем явнее выступают принципиальные особенности, влияющие на возможность его реализации.
На рисунке 1 приведена карта значений периода колебаний в модели Ресслера. Сравнение ее с аналогичной картой для генератора Анищенко-Астахова (рисунок 3 в предыдущей статье) показывает качественное сходство между рассматриваемыми моделями.
На рисунках 2, 3 и 4 приведены характеристики режима для модели Ресслера при a=b=0.2 и c=2. Временные реализации на рисунке 2 показывают, как меняются значения переменных x, y и z с течением времени. Стоит здесь отметить, что здесь время является безразмерным, относительным, введенным специально для этой модели. На рисунке 3 представлена проекция фазового портрета на плоскость переменных x,y. Фазовым портретом называют изображения фазовых траекторий, демонстрирующих поведение модели, в фазовом пространстве. Фазовым пространством называют пространство независимых переменных, в данном случае это пространство трехмерно (переменные x, y и z). Стоит здесь отметить, что с привычным пространством, в котором мы живем это фазовое пространство не имеет ничего общего, фазовое пространство абстрактно и существует только для описания модели, для которой сконструировано. Фазовой траеторией называют множество фазовых состояний, которые получаются друг из друга по некоторому закону эволюции, системой дифференциальных уравнений Ресслера в данном случае. Фазовым состоянием называют любую точку в фазовом пространстве. Если задать любую точку в фазовом пространстве как начальную и применить закон эволюции, из нее стартует фазовая траектория. На рисунке 4 приведен спектр для полученного режима динамики (рисунки 2 и 3), на спектре присутствует гармоника на частоте, обратной периоду колебаний, который нетрудно определеть по временным реализациям на рисунке 2.
Увеличение параметра c приведет к реализации сценария перехода к хаосу через каскад бифуркаций удвоения периода, также, как показано в предыдущей статье.
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
Cкрипты gnuplot для построения картинок
reset
set title "Spectrum for a=b=0.2 and c=2.0"
set rmargin screen 0.99
set lmargin screen 0.13
set bmargin screen 0.15
set tmargin screen 0.85
set xlabel "f" offset 0,0.8
set ylabel "A" rotate by 90 offset 0.6,0
set terminal pngcairo background rgb "white" font "Helvetica,30" fontscale 1.0 size 1000, 600
set output "spectra-0-2-0-2-2-0.png"
set xrange [0.0005:0.18]
set xtics 0.05
set mxtics 5
set ytics ("10^{-5}" 0.00001, "10^{-4}" 0.0001, "10^{-3}" 0.001, "10^{-2}" 0.01,"10^{-1}" 0.1, "10^{0}" 1, "10^{1}" 10)
set log y
plot "spectra-0-2-0-2-2-0.dat" u 3:6 w lp lw 2 pt 7 ps 0.7 notitle
reset
set title "Time realization for a=b=0.2 and c=2.0"
set rmargin screen 0.99
set lmargin screen 0.13
set bmargin screen 0.15
set tmargin screen 0.85
set xlabel "t" offset 0,0.8
set ylabel "x,y,z" rotate by 90 offset 0.6,0
set terminal pngcairo background rgb "white" font "Helvetica,30" fontscale 1.0 size 1000, 600
set output "TS-0-2-0-2-2-0.png"
set key horizontal
set xrange [0.0:49.99]
set xtics 10
set mxtics 5
set yrange [-5:6]
plot "TS-0-2-0-2-2-0.dat" u ($3-10000.0):4 w l lw 2 title "x", "" u ($3-10000.0):5 w l lw 2 title "y", "" u ($3-10000.0):6 w l lw 2 title "z"
reset
set title "Phase portrait for a=b=0.2 and c=2.0"
set rmargin screen 0.99
set lmargin screen 0.13
set bmargin screen 0.15
set tmargin screen 0.85
set xlabel "y" offset 0,0.8
set ylabel "x" rotate by 90 offset 0.6,0
set terminal pngcairo background rgb "white" font "Helvetica,30" fontscale 1.0 size 1000, 600
set output "PP-0-2-0-2-2-0.png"
set xrange [-5:5]
set xtics 2
set mxtics 5
set yrange [-5:5]
set mytics 5
set ytics 2
plot "TS-0-2-0-2-2-0.dat" u 4:5 w l lw 2 notitle
reset
Программа на Си/C++ для вычисления спектров колебаний в модели Ресслера
#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
#include <string>
#include <sstream>
#include <cstdlib>
#include <cstdio>
#include <ctime>
/*
/*The program written by Andrei Bukh, to compile it run:
* $ g++ -o Roessler Roessler.cpp -lm
* $ ./Roessler 0.17 0.17 4.0 0.1 0.1 0.1 20000 30000 10000 0.1
*/
int totalvars = 3;
int totalpars = 3;
int func (double *** vars, int marker, int srfr, int stto, double * times, double * pars) {
double x = vars[marker][srfr][0];
double y = vars[marker][srfr][1];
double z = vars[marker][srfr][2];
double a = pars[0];
double b = pars[1];
double c = pars[2];
vars[marker][stto][0] = -y-z;
vars[marker][stto][1] = x+a*y;
vars[marker][stto][2] = b+z*(x-c);
return 0;
}
int rk4step (double *** vars, double *times, double *pars, int marker, int maxdelay, double h) {
int pm = marker-1; if (pm<0) pm+=maxdelay;
for (int v=0;v<totalvars;v++)
vars[marker][5][v] = vars[marker][0][v] = vars[pm][5][v];
times[marker] = times[pm]; func (vars,marker,0,4,times,pars); /*K1*/
for (int v=0;v<totalvars;v++)
vars[marker][5][v] += h * vars[marker][4][v] / 6.0;
for (int v=0;v<totalvars;v++)
vars[marker][1][v] = vars[marker][0][v] + h * 0.5 * vars[marker][4][v];
times[marker] += 0.5*h; func (vars,marker,1,4,times,pars); /*K2*/
for (int v=0;v<totalvars;v++)
vars[marker][5][v] += h * vars[marker][4][v] / 3.0;
for (int v=0;v<totalvars;v++)
vars[marker][2][v] = vars[marker][0][v] + h * 0.5 * vars[marker][4][v];
func (vars,marker,2,4,times,pars); /*K3*/
for (int v=0;v<totalvars;v++)
vars[marker][5][v] += h * vars[marker][4][v] / 3.0;
for (int v=0;v<totalvars;v++)
vars[marker][3][v] = vars[marker][0][v] + h * vars[marker][4][v];
times[marker] += 0.5*h; func (vars,marker,3,4,times,pars); /*K4*/
for (int v=0;v<totalvars;v++)
vars[marker][5][v] += h * vars[marker][4][v] / 6.0;
return 0;
}
int singleparset (double *trvar, double *trpar, int maxdelay, double time, double trans, double timestep) {
double*** vars; double* times; double* pars;
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];
times[0] = 0.0;
for (int d=0;d<maxdelay;d++)
for (int s=0;s<6;s++)
for (int v=0;v<totalvars;v++) vars[d][s][v] = trvar[v];
for (int p=0;p<totalpars;p++) pars[p] = trpar[p];
std::ofstream ts_file("TS.dat");
int spectra_num = maxdelay/2; if (spectra_num<1) spectra_num = 1;
int spectra_count = 0;
double *spectra_part_im; spectra_part_im = new double[spectra_num];
double *spectra_part_re; spectra_part_re = new double[spectra_num];
double *spectra_im; spectra_im = new double[spectra_num];
double *spectra_re; spectra_re = new double[spectra_num];
int marker = 0;
long int t = 0;
double per_x_test, per_y_test, per_z_test; double per_start_x, per_start_y, per_start_z; double per_start_time=trans; double per_finish_time=time; bool per_started = false; bool per_gone = false; bool per_near_target = false; bool per_finished = false; double mindist = 0.01;
for (double ct=0.0;ct<time;ct+=timestep) {
t++;
marker++; if (marker>=maxdelay) marker = 0;
int prevmarker = marker-1; if (prevmarker<0) prevmarker = maxdelay-1;
rk4step (vars, times, pars, marker, maxdelay, timestep);
if (((t%(1000000))==0)) std::cout<<"ct="<<ct<<"\n";
if (ct>trans) {
if (!per_started) {
per_started = true;
per_start_time = ct;
per_start_x = vars[marker][0][0];
per_start_y = vars[marker][0][1];
per_start_z = vars[marker][0][2];
}
else {
if (!per_gone) {
if ((ct-per_start_time)>(3.0)) {
per_gone = true;
}
}
else {
if (!per_finished) {
double dx = (vars[marker][0][0]-per_start_x);
double dy = (vars[marker][0][1]-per_start_y);
double dz = (vars[marker][0][2]-per_start_z);
double dist_loc = sqrt(dx*dx+dy*dy+dz*dz);
if (!per_near_target) {
if (dist_loc<mindist) {
mindist = dist_loc;
per_finish_time = ct;
per_near_target = true;
}
}
else {
if (dist_loc<mindist) {
mindist = dist_loc;
per_finish_time = ct;
}
else {
per_finished = true;
}
}
}
}
}
if (((t%(10000))==0)) std::cout<<"ct="<<ct<<"\n";
ts_file << pars[0] << " " << pars[1] << " " << ct << " " << vars[marker][0][0] << " " << vars[marker][0][1] << " " << vars[marker][0][2] << "\n";
if (marker==(maxdelay-1)) {
if (spectra_num>2) {
for (int k=0; k<spectra_num; k++) {
spectra_part_re[k] = 0.0;
spectra_part_im[k] = 0.0;
for (int j=0;j<maxdelay;j++) {
spectra_part_re[k] += vars[j][0][0] * cos(2.0 * M_PI * k * j / maxdelay);
spectra_part_im[k] -= vars[j][0][0] * sin(2.0 * M_PI * k * j / maxdelay);
}
spectra_re[k]+=spectra_part_re[k]/maxdelay;
spectra_im[k]+=spectra_part_im[k]/maxdelay;
}
spectra_count++;
}
}
}
}
ts_file.close();
double maxamp = 0.0; double freqmax = 0.0;
if (spectra_num>2) {
std::ofstream spec("spectra.dat");
for (int k=0; k<spectra_num; k++) {
double amp = sqrt(spectra_re[k]*spectra_re[k]+spectra_im[k]*spectra_im[k]);
double freq = 1.0*k/(maxdelay*timestep);
if (amp>maxamp) {
freqmax = freq;
maxamp = amp;
}
spectra_re[k]/=spectra_count;
spectra_im[k]/=spectra_count;
spec << pars[0] << " " << pars[1] << " " << freq << " " << spectra_re[k] << " " << spectra_im[k] << " " << amp << " " << atan(spectra_im[k]/spectra_re[k]) << "\n";
}
spec.close();
}
std::ofstream specmax("spectra_max.dat",std::ios::app);
specmax << pars[0] << " " << pars[1] << " " << pars[2] << " " << freqmax << " " << maxamp << " " << per_finish_time-per_start_time << " " << mindist << " " << std::setprecision(25) << vars[marker][0][0] << " " << std::setprecision(25) << vars[marker][0][1] << " " << std::setprecision(25) << vars[marker][0][2] << "\n";
specmax.close();
for (int v=0;v<totalvars;v++) trvar[v] = vars[marker][0][v];
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;
}
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;
}
return res;
}
int main (int argc, char* argv[]) {
double pars[totalpars];
pars[0] = takefloatpar(argc, argv, 1);
pars[1] = takefloatpar(argc, argv, 2);
pars[2] = takefloatpar(argc, argv, 3);
double vars[totalvars];
vars[0] = takefloatpar(argc, argv, 4);
vars[1] = takefloatpar(argc, argv, 5);
vars[2] = takefloatpar(argc, argv, 6);
int maxdelay = takeintpar(argc, argv, 7);
double time = takefloatpar(argc, argv, 8);
double trans = takefloatpar(argc, argv, 9);
double timestep = takefloatpar(argc, argv, 10);
std::ifstream ICf ("IC.dat");
if (ICf.is_open()) {
if (!ICf.eof()) ICf >> vars[0];
if (!ICf.eof()) ICf >> vars[1];
if (!ICf.eof()) ICf >> vars[2];
ICf.close();
}
singleparset(vars,pars,maxdelay,time,trans,timestep);
return 0;
}
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
