LoginSignup
4
2

More than 3 years have passed since last update.

@Amefurase(Ame furase)

Processingでカオスを感じよう

Last updated at Posted at 2019-09-15

はじめに

カオスとは何か？

カオスについて詳しくないので、自分なりに「初期値のごくわずかなずれが指数関数的に増大していくために、有限な桁の計算しか行えないコンピュータでは長期的な予測が不可能な力学系」と考えています。生活中での例でいうと気象予報などが、それにあたるため長期的な予測が不可能なようです。[ローレンツ方程式]

使用する言語

映像の描画が簡単な言語で、友人が使っていたのでProcessingを選択しました。

今回対象にするカオス

カオスを説明する際によく例に出され、挙動が面白い二重振り子を選択しました。

準備

  • ここからprocessingをインストール

  • IDEのフォント設定を変更し、コメントアウトを有効化[参考URL]

二重振り子の軌道計算の内容

運動方程式はこちらの記事での座標系の運動方程式を使いました。
微分方程式の数値計算には、運動方程式を1階化して4次のルンゲクッタ法を使いました。

アニメーション

1個の振子のアニメーション

コード

コードの大まかな流れとしては、以下のようになります。
1. 二重振り子クラスdpを定義
2. dpクラスのオブジェクトfurikoを宣言
3. setup関数やsettings関数でウィンドウサイズ、フレームレートを設定
4. 各フレーム毎に実行されるdraw関数の中で、furikoオブジェクトの計算および描画関数を実行
5. 以降4を繰り返す

double_pendulum.pde
//ウィンドウサイズ
int window = 800;
//フレームレート
int fps = 30;
//円周率
float pi = 3.141592;

//二重振り子オブジェクトを定義
dp furiko;

void settings() {
  size(window,window);
}
void setup(){
  background(0);
  frameRate(fps);
  //二重振り子のコンストラクタ
  furiko = new dp(pi *0.95 ,pi * 0.99,window / 2,window / 2);
}

void draw(){
  background(0);
  //二重振り子のリンク角度を計算
  furiko.draw();
  saveFrame("frames/####.png");
}

//2重振り子クラスを定義
class dp {
  //支点の座標
  int x,y;
  //第一角度、第二角度、泰一角速度、第二角速度、刻み幅
  float th1,th2,dth1,dth2,ddth1,ddth2,dt;
  //第一振子長さ、第二振子長さ,重り質量1,2,重力加速度
  float r1 = 2;
  float r2 = 1;
  float m1 = 1;
  float m2 = 1;
  float g = 9.81;

  //ルンゲクッタの係数
  float[] kt1 = new float[4];
  float[] kt2 = new float[4];
  float[] kd1 = new float[4];
  float[] kd2 = new float[4];

  //コンストラクタ
  dp(float th1_,float th2_,int x_,int y_){
    th1 = th1_;
    th2 = th2_;
    dth1 = 0.0;
    dth2 = 0.0;
    //fpsから計算の刻み幅決定
    dt = 1.0 / fps;
    x = x_;
    y = y_;
  }
  //リンク角度を計算
  void calc(){
    //第1次近似を計算
    kt1[0] = dt * Fth1(th1,th2,dth1,dth2);
    kt2[0] = dt * Fth2(th1,th2,dth1,dth2);
    kd1[0] = dt * Fdth1(th1,th2,dth1,dth2);
    kd2[0] = dt * Fdth2(th1,th2,dth1,dth2);
    //第2次近似以降を計算
    for(int i = 1;i <= 3; i++){
      if(i != 3){
        kt1[i] = dt * Fth1(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
        kt2[i] = dt * Fth2(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
        kd1[i] = dt * Fdth1(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
        kd2[i] = dt * Fdth2(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
      }else{
        kt1[i] = dt * Fth1(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
        kt2[i] = dt * Fth2(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
        kd1[i] = dt * Fdth1(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
        kd2[i] = dt * Fdth2(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
      }
    }
    //dt後の座標を更新
    th1 += (kt1[0] + 2 * kt1[1] + 2 * kt1[2] + kt1[3]) /6.0;
    th2 += (kt2[0] + 2 * kt2[1] + 2 * kt2[2] + kt2[3]) /6.0;
    dth1 += (kd1[0] + 2 * kd1[1] + 2 * kd1[2] + kd1[3]) /6.0;
    dth2 += (kd2[0] + 2 * kd2[1] + 2 * kd2[2] + kd2[3]) /6.0;
  }
  //第一角速度を計算
  float Fth1(float t1,float t2,float d1,float d2){
      return d1;
  }
  //第二角速度を計算
  float Fth2(float t1,float t2,float d1,float d2){
      return d2;
  }
  //第一角加速度を計算
  float Fdth1(float t1,float t2,float d1,float d2){
    float dd1;
    dd1 = (g * m1 * sin(t1) + g * m2 * sin(t1 - 2 * th2) / 2.0 +g * m2 *sin(t1) / 2.0 + m2 * r1 * sin(2 *t1- 2 *t2) * pow(d1,2) / 2.0  + m2 * r2 * sin(t1 - t2) * pow(d2 ,2)) / (r1 *(-m1 + m2 * pow(cos(t1 - t2) , 2) -m2));
    return dd1;
  }
  //第二角加速祖を計算
  float Fdth2(float t1,float t2,float d1,float d2){
    float dd2;
    dd2 = ((m1 + m2) * (g * sin(t2) - r1 * sin(t1 - t2) *pow(d1,2) ) - (g * m1 * sin(t1) + g * m2 *sin(t1) + m2 * r2 * sin(t1 - t2) * pow(d2 , 2)) * cos(t1 - t2)) / (r2 * (-m1 + m2 * pow(cos(t1 - t2),2) - m2));
    return dd2;
  }
  //描画関数
  void draw(){
    //
    calc();
    float x1,y1,x2,y2;
    //リンク端の座標を計算
    x1 = x + r1 * window/8 * cos(-th1 + pi/2);
    y1 = y + r1 * window/8 * sin(-th1 + pi/2);
    x2 = x1 + r2 * window/8 * cos(-th2 + pi/2);
    y2 = y1 + r2 * window/8 * sin(-th2 + pi/2);
    //リンクを描画
    stroke(255);
    line(x,y,x1,y1);
    line(x1,y1,x2,y2);
    //リンク端を描画
    fill(255);
    ellipse(x1,y1,5,5);
    ellipse(x2,y2,5,5);
  }
}
}

アニメーション

furiko.gif

二重振り子のアニメーションが描画できました！
しかしながらこれだけでは、この運動がカオス的な振る舞いをしているかわからないので、ここから一気に振り子を100個に拡張していきます。難しそうに感じますが、2重振り子をクラスで定義しているので変数を変えるだけでOKです。（始めからクラスで定義しなくてもいいですが、実装出来たらクラスで定義し直しておく拡張性があがるので便利です。）

100個の振子のアニメーション

$10^{-5}$[rad]だけずらした、100個の二重振り子のアニメーション

コード

dpクラスのオブジェクトを配列で生成する。

n_double_pendulum.pde
//ウィンドウサイズ
int window = 800;
//フレームレート
int fps = 30;
//円周率
float pi = 3.141592;
//振子の数
int n =100;
//振子のずらす角度
float dtheat = 0.00001;

//振子オブジェクトを定義
dp furiko[] = new dp[n];

void settings() {
  size(window,window);
}

void setup(){
  background(0);
  frameRate(fps);
  //振子のコンストラクタ
 float theat = 0;
  for(int i = 0;i<n;i++){
    furiko[i] = new dp(pi *0.95 ,pi * 0.99 +theat,window / 2,window / 2);
    theat += dtheat;
  }
}

void draw(){
  background(0);
  //リンク角度を計算
  for(int i = 0;i<n;i++){
    furiko[i].draw();
  }
  //コメントアウトを外すと角フレームを保存して、アニメーション化
  saveFrame("frames30fps/####.png");
}

//2重振り子クラスを定義
class dp {
  //支点の座標
  int x,y;

  //第一角度、第二角度、泰一角速度、第二角速度、刻み幅
  float th1,th2,dth1,dth2,ddth1,ddth2,dt;
  //第一振子長さ、第二振子長さ,重り質量1,2,重力加速度
  float r1 = 2;
  float r2 = 1;
  float m1 = 0.5;
  float m2 = 0.5;
  float g = 9.81;
  //ルンゲクッタの係数
  float[] kt1 = new float[4];
  float[] kt2 = new float[4];
  float[] kd1 = new float[4];
  float[] kd2 = new float[4];
  //コンストラクタ
  dp(float th1_,float th2_,int x_,int y_){
    th1 = th1_;
    th2 = th2_;
    dth1 = 0.0;
    dth2 = 0.0;
    //fpsから計算の刻み幅決定
    dt = 1.0 / fps;
    x = x_;
    y = y_;
  }

  //リンク角度を計算
  void calc(){
    //第1次近似を計算
    kt1[0] = dt * Fth1(th1,th2,dth1,dth2);
    kt2[0] = dt * Fth2(th1,th2,dth1,dth2);
    kd1[0] = dt * Fdth1(th1,th2,dth1,dth2);
    kd2[0] = dt * Fdth2(th1,th2,dth1,dth2);
    //第2次近似以降を計算
    for(int i = 1;i <= 3; i++){
      if(i != 3){
        kt1[i] = dt * Fth1(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
        kt2[i] = dt * Fth2(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
        kd1[i] = dt * Fdth1(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
        kd2[i] = dt * Fdth2(th1 + kt1[i-1] / 2,th2 + kt2[i-1] /2,dth1 + kd1[i-1] /2,dth2 +kd2[i-1] /2);
      }else{
        kt1[i] = dt * Fth1(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
        kt2[i] = dt * Fth2(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
        kd1[i] = dt * Fdth1(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
        kd2[i] = dt * Fdth2(th1 + kt1[i-1],th2 + kt2[i-1] ,dth1 + kd1[i-1] ,dth2 +kd2[i-1] );
      }
    }
    //dt後の座標を更新
    th1 += (kt1[0] + 2 * kt1[1] + 2 * kt1[2] + kt1[3]) /6.0;
    th2 += (kt2[0] + 2 * kt2[1] + 2 * kt2[2] + kt2[3]) /6.0;
    dth1 += (kd1[0] + 2 * kd1[1] + 2 * kd1[2] + kd1[3]) /6.0;
    dth2 += (kd2[0] + 2 * kd2[1] + 2 * kd2[2] + kd2[3]) /6.0;
  }

  //第一角速度を計算
  float Fth1(float t1,float t2,float d1,float d2){
      return d1;
  }
  //第二角速度を計算
  float Fth2(float t1,float t2,float d1,float d2){
      return d2;
  }
  //第一角加速度を計算
  float Fdth1(float t1,float t2,float d1,float d2){
    float dd1;
    dd1 = (g * m1 * sin(t1) + g * m2 * sin(t1 - 2 * th2) / 2.0 +g * m2 *sin(t1) / 2.0 + m2 * r1 * sin(2 *t1- 2 *t2) * pow(d1,2) / 2.0  + m2 * r2 * sin(t1 - t2) * pow(d2 ,2)) / (r1 *(-m1 + m2 * pow(cos(t1 - t2) , 2) -m2));
    return dd1;
  }
  //第二角加速祖を計算
  float Fdth2(float t1,float t2,float d1,float d2){
    float dd2;
    dd2 = ((m1 + m2) * (g * sin(t2) - r1 * sin(t1 - t2) *pow(d1,2) ) - (g * m1 * sin(t1) + g * m2 *sin(t1) + m2 * r2 * sin(t1 - t2) * pow(d2 , 2)) * cos(t1 - t2)) / (r2 * (-m1 + m2 * pow(cos(t1 - t2),2) - m2));
    return dd2;
  }

  //描画関数
  void draw(){
    calc();
    float x1,y1,x2,y2;
    //リンク端の座標を計算
    x1 = x + r1 * window/8 * cos(-th1 + pi/2);
    y1 = y + r1 * window/8 * sin(-th1 + pi/2);
    x2 = x1 + r2 * window/8 * cos(-th2 + pi/2);
    y2 = y1 + r2 * window/8 * sin(-th2 + pi/2);
    //リンクを描画
    stroke(255);
    line(x,y,x1,y1);
    line(x1,y1,x2,y2);
    //リンク端を描画
    fill(255);
    ellipse(x1,y1,5,5);
    ellipse(x2,y2,5,5);
  }
}

アニメーション

nfuriko30fps2.gif

初期には変化のなかった振子が突然、振り子がばらばらに動き始めるます。
このことから初期値に微小な変化があったときに、長期的な予測が不可能になるというカオス的振る舞いを確認できました。

おわりに

Processingは映像の描画や環境を用意するのが簡単なので、いろいろなアニメーションに挑戦していきたいと思います。
（追記）
振り子が1個の時と100個の時で質量のパラメータが違うため、運動の速度が異なっています。

4
2
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
Sign upLogin
4
2