二重振り子はカオスな運動をすることで知られている。その軌跡をProcessingで書くと美しい。
単振り子の作り方
The Nature of Code「3.9 Trigonometry and Forces: The Pendulum」を読もう。
http://natureofcode.com/book/chapter-3-oscillation/
二重振り子の加速度
\theta_1'' = \frac{-g(2m_1 + m_2){\sin \theta_1} - m_2g{\sin (\theta_1 - 2\theta_2) -2{\sin (\theta_1 - \theta_2)}}m_2(\theta_2'^2L_2+\theta_1'^2L_1{\cos(\theta_1 - \theta_2))}}{L_1(2m_1+m_2-m_2{\cos(2\theta_1 - 2\theta_2))}}
\theta_2'' = \frac{2{\sin (\theta_1 - \theta_2)(\theta_1'^2L_1(m_1+m_2)+g(m_1+m_2){\cos \theta_1}}+\theta_2'^2L_2m_2{\cos(\theta_1 - \theta_2))}}{L_2(2m_1+m_2-m_2{\cos(2\theta_1 - 2\theta_2))}}
二重振り子の加速度は上のようなエグい数式で表される。
親の振り子
内側の振り子を親とする。
角速度
角度の一回微分が角速度になることに注意する。
float alpha = -gravity * (2 * mass1 + mass2) * sin(angle)-mass2 * gravity * sin(angle - 2 * childAngle) - 2 * sin(angle - childAngle) * mass2 *(childAngleVelocity * childAngleVelocity * l2 + angleVelocity * angleVelocity * l1 * cos(angle - childAngle));
float beta = l1 * (2 * mass1 + mass2 - mass2 * cos(2 * angle - 2 * childAngle));
angleAcceleration = alpha / beta;
角度を求める
angleVelocity += angleAcceleration;
angle += angleVelocity;
子の振り子
外側の振り子を子とする。
角速度
float gunnma = 2 * sin(angle - childAngle) * (angleVelocity * angleVelocity * l1 *(mass1 + mass2) + gravity * (mass1 + mass2) * cos(angle) + childAngleVelocity * childAngleVelocity * l2 * mass2 * cos(angle - childAngle));
float phai = l2 * (2 * mass1 + mass2 - mass2 * cos(2 * angle - 2 * childAngle));
childAngleAcceleration = gunnma / phai;
角度を求める
childAngleVelocity += childAngleAcceleration;
childAngle += childAngleVelocity;
描画
void display() {
// 親の振り子
location.set(l1 * sin(angle), l1 * cos(angle), 0);
location.add(origin);
translate(0, height / 3);
float c = map(location.x,0, width, 0, 255);
stroke(96, 230,c);
fill(255);
line(origin.x, origin.y, location.x, location.y);
ellipse(location.x, location.y, mass1 * 3, mass1 * 3);
// 子の振り子
childLocation = location.get();
childLocation.add(l2 * sin(childAngle), l2 * cos(childAngle));
line(location.x, location.y, childLocation.x, childLocation.y);
ellipse(childLocation.x, childLocation.y, mass2 * 3, mass2 * 3);
// ellipse(childLocation.x, childLocation.y, 1, 1);
}