はじめに
離散写像の1次元分岐図は簡単なプログラムで実装できる.
これまでは必要に応じて昔のプログラムを参照していたが,ファイルを探すのが面倒であったりするので Qiita 上に残すことにした.
必要に応じて適宜修正を加えていく.
現状,汎用性の低い書き方になっており,あまりおもしろくないので,以降の修正で検討していく.
Logistic map
Logistic map は,代表的な非線形力学系であり,次式で定義される.
$$
x_{n+1} = rx_n(1-x_n)
$$
力学系の解軌道 $\{x_0, x_1, x_2, \ldots\}$ はあるパラメタを境にカオス的に振る舞う.
C 言語によるコード
ここでは C 言語で実装したが,どの言語で書くにしてもプログラムの大枠は変わらない.
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
// Definition of Logistic map
double logistic_map(double x, double r){
return (r*x*(1.0-x));
}
int main (void){
double range[2] = {2.4, 4.0}; // Parameter domain
double r = range[0], x0 = 0.5; // Initial parameter and state
int resolution = 2000, maps= 2000; // Resolution of diagram, count of maps
double step = (range[1]-range[0])/(double)resolution;
FILE *fp;
fp = fopen("out.bf1", "w"); // Output filename
// Print initial value, parameter, and environment information
char header[256];
sprintf(header,
"# Logistic map bifurcation diagram\n"
"# x0 = %lf, range = [%lf %lf], resolution = %d, maps = %d\n",
x0, range[0], range[1], resolution, maps);
fprintf(fp, "%s", header);
fprintf(stdout, "%s", header);
// Main loop
for(int ri = 0; ri <= resolution; ri++){
fprintf(fp, "%lf %lf\n", r, x0);
for(int mi = 1; mi <= maps; mi++){
x0 = logistic_map(x0, r);
// if x0 got divergent
if (isinf(x0)){
fprintf(stderr, "Got divergent at ri = %d, mi = %d\n", ri, mi);
fclose(fp);
return -1;
}
// else print
fprintf(fp, "%lf %lf\n", r, x0);
fprintf(stdout, "\t[%04d/%04d] [%04d/%04d] %+.5lf %+.5lf\r",
ri, resolution, mi, maps, r, x0);
}
// iterate paramter
r += step;
fprintf(fp, "\n");
}
fprintf(stdout, "\n");
fclose(fp);
return 0;
}
Gnuplot による図示
上記のプログラムによって出力される out.bf1
から,次のような図が得られる.
plot のいろはは別の記事に任せることにする.