processingの備忘録　-二重振り子- - プログラミングの備忘録

こんにちは。
今回は「二重振り子」をprocessingで描画してみたいと思います。

二重振り子をprocessingでやるのはn番煎じ的なところがありますが、やりたいのでやります。

二重振り子とは
プログラム化する
おまけ
まとめ
参考

二重振り子とは

「二重振り子」はその名の通り、振り子が2つつながったものです。
カオスの例としてよく挙げられるので、聞いたり見たりしたことがあるかもしれません。

二重振り子において、各振り子の角度の二階微分については以下のように表されるようです。

$$ \begin{align} \ddot{\theta}_1 &= \frac{ -m_2 ~ l_2 ~ \ddot{\theta}_2 ~ \mathrm{cos}(\theta_1 - \theta_2) - m_2 ~ l_2 ~ \dot{\theta}_2^2 ~ \mathrm{sin}(\theta_1 - \theta_2) - (m_1+m_2) ~ g ~ \mathrm{sin}\theta_1 }{ (m_1+m_2) ~ l_1 } \\ \ddot{\theta}_2 &= \frac{ m_2 ~ l_1 ~ \dot{\theta}_1^2 ~ \mathrm{sin}(\theta_1 - \theta_2) - m_2 ~ l_1 ~ \ddot{\theta}_1 ~ \mathrm{cos}(\theta_1 - \theta_2) - m_2 ~ g ~ \mathrm{sin}\theta_2 }{ m_2 ~ l_2 } \end{align} $$

この数式の導出については「参考」の項にあるサイトなどを参照してみてください。

プログラム化する

さっそくつくっていきます。

先程の式で「角度の変化量の変化量ddtheta」が求められるので、これに時間の変化量をかければ「角度の変化量dtheta」になり、さらに時間の変化量をかければ「角度theta」になります。

ここから振り子の重り部分の座標 $(x, y)$ が求められるので、これを描画することで二重振り子が完成する、といった感じです。

float ddtheta1, ddtheta2;
float dtheta1, dtheta2;
float theta1, theta2;
float m1, m2;
float l1, l2;
float x1, y1, x2, y2;
float g;
float dt;

void setup(){
  size(500, 500);
  
  //振り子の角度
  theta1 = PI/2;
  theta2 = PI/2;
  
  //重りの重さ
  m1 = 10;
  m2 = 1;
  
  //腕の長さ
  l1 = 100;
  l2 = 100;
  
  //重力加速度
  g = 10;
  
  //時間の進む速さ
  dt = 0.2;
}

void draw(){
  background(0);
  
  //角度の計算
  ddtheta1 = (- m2*l2*ddtheta2*cos(theta1-theta2) - m2*l2*dtheta2*dtheta2*sin(theta1-theta2) - (m1+m2)*g*sin(theta1))/((m1+m2)*l1);
  ddtheta2 = (m2*l1*dtheta1*dtheta1*sin(theta1-theta2) - m2*l1*ddtheta1*cos(theta1-theta2) - m2*g*sin(theta2))/(m2*l2);
  
  dtheta1 += dt * ddtheta1;
  dtheta2 += dt * ddtheta2;
  
  theta1 += dt * dtheta1;
  theta2 += dt * dtheta2;
  
  //座標の計算
  x1 = l1*sin(theta1);
  y1 = l1*cos(theta1);
  
  x2 = l1*sin(theta1) + l2*sin(theta2);
  y2 = l1*cos(theta1) + l2*cos(theta2);

  //描画
  translate(width/2, height/4);
  noStroke();
  ellipse(x1, y1, 20, 20);
  ellipse(x2, y2, 20, 20);
  stroke(255);
  line(0, 0, x1, y1);
  line(x1, y1, x2, y2); 
}

式がわかれば後は簡単ですね。

実行すると、

よく見る感じです。

おまけ

角度を0.00001 radずつ変えて100個並べたものを動かしてみました。

初めは1つに見えますが、しばらくすると離れてきて最終的にはばらばらになってしまいました。

突貫工事ですが、コードも載せておきます。

float g, dt;
int n = 100;
DoublePendulum[] pendulums = new DoublePendulum[n];

void setup(){
  size(500, 500);
  
  //振り子の角度
  for(int i = 0; i < n; i ++){
    pendulums[i] = new DoublePendulum(PI/3+i*0.00001, PI/2, map(i, 0, n, 100, 255));
  }
  
  g = 10;
  
  dt = 0.2;

}

void draw(){
  background(0);
  translate(width/2, height/3);
  
  for(int i = 0; i < n; i ++){
    pendulums[i].move();
    pendulums[i].display();
  }
}

class DoublePendulum{
  float ddtheta1, ddtheta2;
  float dtheta1, dtheta2;
  float theta1, theta2;
  float m1, m2;
  float l1, l2;
  float x1, y1, x2, y2;
  float c;

  DoublePendulum(float _t1, float _t2, float _c){
    theta1 = _t1;
    theta2 = _t2;
    
    m1 = 10;
    m2 = 1;
    
    l1 = 100;
    l2 = 100;
    
    c = _c; 
  }

  void move(){
    ddtheta1 = (- m2*l2*ddtheta2*cos(theta1-theta2) - m2*l2*dtheta2*dtheta2*sin(theta1-theta2) - (m1+m2)*g*sin(theta1))/((m1+m2)*l1);
    ddtheta2 = (m2*l1*dtheta1*dtheta1*sin(theta1-theta2) - m2*l1*ddtheta1*cos(theta1-theta2) - m2*g*sin(theta2))/(m2*l2);
    
    dtheta1 += dt * ddtheta1;
    dtheta2 += dt * ddtheta2;
    
    theta1 += dt * dtheta1;
    theta2 += dt * dtheta2;
    
    x1 = l1*sin(theta1);
    y1 = l1*cos(theta1);
    
    x2 = l1*sin(theta1) + l2*sin(theta2);
    y2 = l1*cos(theta1) + l2*cos(theta2);
  }

  void display(){
    noStroke();
    fill(c);
    ellipse(x1, y1, 20, 20);
    ellipse(x2, y2, 20, 20);
    stroke(c);
    line(0, 0, x1, y1);
    line(x1, y1, x2, y2); 
  }
}