0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

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

Last updated at Posted at 2021-02-18

はじめに

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

out.jpg

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?