![Рисунок 1. Фазовый портрет генератора Анищенко-Астахова при m=1.563, g=0.17 [Результат авторов статьи https://fizika.sgu.ru/system/files_force/text-pdf/2020/01/2005_1_54_anishchenko.pdf?download=1].](https://avatars.dzeninfra.ru/get-zen_doc/271828/pub_65d3079e83a4b769b7d1aa94_65d32b104f3f3438f6d59fe7/scale_1200)
Типичный вид хаотического аттрактора приведен на рисунке 1. Уравнения показаны на рисунке 2.
Генератор Анищенко-Астахова является важной моделью нелинейной динамике, в которой реализуется возникновение хаотических автоколебаний, в том числе через каскад бифуркаций удвоения периода, который рассматривается в этой статье.
На рисунке 3 приведена карта значений периода колебаний в отдельном генераторе Анищенко-Астахова. Базовым режимом на этой карте является цикл периода 1 (черная область). Если фиксировать один из параметров и увеличивать другой, можно наблюдать каскад удвоения периода. В этой статье обсудим еще одну характеристику (спектр), позволяющую наглядно показать переход от периодических колебаний к хаосу через каскад бифуркаций удвоения периода. Наблюдать бифуркации удвоения периода можно с помощью значений мультипликаторов, величины показателей Ляпунова, временных реализаций, фазовых портретов и спектров. Последние будут отражать появление новых кратных гармоник при удвоении колебаний.
На рисунке 4 приведены спектры колебаний генератора Анищенко-Астахова при фиксированном параметре g=0.17 и увеличивающихся значениях параметра m. Видно, что при m=1 и m=1.08 в области низких частот наблюдается только одна значительная по амплитуде гармоника, близкая к значению f=0.16, соответствующему значению T=6.25. При m=1.09 в спектре возникает дополнительная гармоника на половинной частоте по отношению к основной. Увеличение m до значения m=1.2 и до значения m=1.44 (близкого к следующей бифуркации) качественно не изменяет спектр колебаний в области низких частот. После очередной бифуркации при m=1.45 в рассматриваемом диапазоне частот в спектре появляется еще две гармоники, одна из которых (самая нижняя) соответствует дважды удвоенному периоду. При m=1.48 и при m=1.51 набор гармоник не изменился. Новые гармоники в спектре появляются при m=1.52, индикатором очередного удвоения является опять наиболее низкая гармоника. При m=1.53 набор гармоник почти сохраняется, но самая нижняя остается прежней. Изменение набора гармоник в данном случае связано с низким разрешением спектра, которое можно увеличить, изменив 6-й параметр при вызове программы, приведенной в конце статьи. Наконец, при превышении критического значения параметра m, зависящего от g, генератор переходит в режим хаоса (см. пример спектра при m=1.55 на рисунке 4).
p. s. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
Скрипт для построения спектров в gnuplot
set title "Spectrum for m=1.08 and g=0.17"
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.png"
set xrange [0:0.18]
set xtics 0.05
set mxtics 5
# set yrange [0.0001:20]
set ytics ("10^{-4}" 0.0001, "10^{-3}" 0.001, "10^{-2}" 0.01,"10^{-1}" 0.1, "10^{0}" 1, "10^{1}" 10)
# set mytics 10
# set xtics 0.3
# set ytics 0.1
# set log x
set log y
plot "spectra.dat" u 3:6 w lp lw 2 pt 7 ps 0.7 notitle
Программа для построения спектров по временным реализациям генератора Анищенко-Астахова
#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 GIN GIN.cpp -lm
* $ ./GIN 1.563 0.17 0.1 0.1 0.1 20000 30000 10000 0.1
*/
int totalvars = 3;
int totalpars = 2;
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 m = pars[0];
double g = pars[1];
vars[marker][stto][0] = m*x+y-x*z;
vars[marker][stto][1] = -x;
vars[marker][stto][2] = -g*z+0.5*g*x*(x+fabs(x));
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;
}
}
}
// else {
// 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 (dist_loc<(mindist*0.01)) {
// mindist = dist_loc;
// per_finish_time = ct;
// }
// }
}
}
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] << " " << 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);
double vars[totalvars];
vars[0] = takefloatpar(argc, argv, 3);
vars[1] = takefloatpar(argc, argv, 4);
vars[2] = takefloatpar(argc, argv, 5);
int maxdelay = takeintpar(argc, argv, 6);
double time = takefloatpar(argc, argv, 7);
double trans = takefloatpar(argc, argv, 8);
double timestep = takefloatpar(argc, argv, 9);
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. Чтобы сразу увидеть новый материал в моем блоге в своей ленте, подписывайтесь! Буду рад комментариям, вопросам, предложениям.
