C 言語による Logistic map の1次元分岐図計算プログラム

10281 ワード
C C テキストリンク

はじめに

離散写像の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 のいろはは別の記事に任せることにする．