PendulumTrace ができるまで — 二重振り子を RK4 で解く

PendulumTrace は、二重振り子のカオスな軌跡を描くシミュレーターです。ラグランジュの運動方程式を一階の連立常微分方程式に直し、4 次ルンゲ・クッタ法で積分しています。計算が壊れていないかを全エネルギーで検算しながら動かしているのが、このサービスで一番こだわった部分です。

#🎯 二重振り子が題材として面白い理由

振り子を 2 つ繋いだだけなのに、その動きは解析的に解けません。初期条件をほんのわずか変えるだけで、数秒後にはまったく違う軌跡になります。カオスの教科書的な例です。

そして実装の観点でも都合が良い題材でした。

• 状態が 4 つの数(角度 2 つ、角速度 2 つ)だけで表せる
• なのに出てくる絵が毎回違う
• 正しく計算できているかを検算できる(後述)

#📐 運動方程式を一階に落とす

ラグランジアンから出る二重振り子の運動方程式は、角度の二階微分方程式です。数値積分しやすいよう、一階の連立に書き換えます。

やることは単純で、角速度 ω を状態に加えるだけです。

dθ1/dt = ω1
dθ2/dt = ω2
dω1/dt = (角度と角速度から決まる式)
dω2/dt = (同上)

これで「4 つの変数の、それぞれの時間微分」という形になり、あとは常微分方程式の解法がそのまま使えます。

実際の加速度の式はそこそこ長いです。

const den = 2 * m1 + m2 - m2 * Math.cos(2 * d)
const a1 =
  (-g * (2 * m1 + m2) * Math.sin(theta1) -
    m2 * g * Math.sin(theta1 - 2 * theta2) -
    2 * sd * m2 * (omega2 * omega2 * L2 + omega1 * omega1 * L1 * cd)) /
  (L1 * den)

長いですが、書き写すだけです。ここを自分で導出しようとすると事故ります。既知の式を正確に写して、正しさは後述の検算で担保する、という方針にしました。

#🔁 4 次ルンゲ・クッタ

積分方法は RK4 です。オイラー法だと二重振り子はすぐに破綻します(エネルギーがどんどん増えて、振り子が回り続ける)。

RK4 は 1 ステップの中で微分を 4 回評価して、重み付き平均を取ります。

const k1 = derivatives(s, p)
const k2 = derivatives(scale(k1, dt / 2, s), p)
const k3 = derivatives(scale(k2, dt / 2, s), p)
const k4 = derivatives(scale(k3, dt, s), p)
return {
  theta1: s.theta1 + (dt / 6) * (k1.theta1 + 2*k2.theta1 + 2*k3.theta1 + k4.theta1),
  ...
}

始点・中点 2 回・終点で傾きを測って、(k1 + 2k2 + 2k3 + k4) / 6 で混ぜる。中点を 2 倍に重み付けするのがシンプソン則と同じ発想です。

実効 240Hz(1 秒を 240 ステップ)で回しています。60fps の描画に対して 1 フレームあたり 4 ステップ。刻みを粗くすると、速く振れている瞬間に軌跡が角ばります。

#✅ 全エネルギーで検算する

ここが一番大事だと思っている部分です。数値積分は「それらしく動いているのに間違っている」ことがあります。

二重振り子は(減衰なしなら)全エネルギーが保存します。運動エネルギーと位置エネルギーの和が一定になるはずです。

export function totalEnergy(s: PendulumState, p: PendulumParams): number {
  const T =
    0.5 * m1 * L1 * L1 * omega1 * omega1 +
    0.5 * m2 * (L1*L1*omega1*omega1 + L2*L2*omega2*omega2 +
      2 * L1 * L2 * omega1 * omega2 * Math.cos(theta1 - theta2))
  const V = ...
}

これを画面に出しています。動かしている間、この数字がほぼ一定なら計算は合っている。 じわじわ増減したら、刻みが粗すぎるか式を写し間違えているかのどちらかです。

実際、開発中に運動エネルギーの交差項(2 * L1 * L2 * ω1 * ω2 * cos(θ1-θ2))を落としていて、それはエネルギー表示が暴れることで気づきました。答え合わせの手段を先に用意しておくと、間違いが自分から出てきます。

これは他のサービスでも同じ考え方をしていて、SQL道場では全演習の模範解答を実際に実行して判定を通るか確かめています。「数字を出すものは、検算できる形にしてから作る」という線引きです。

#🎨 軌跡の描き方

下端の球の位置を毎フレーム記録して、線で繋いでいきます。ポイントは 2 つ。

• 古い軌跡を薄くしていく。全部同じ濃さだと、すぐ画面が埋まります
• 座標変換は最後に一度だけ。物理は原点をピボットにした実座標で計算し、描画のときにだけ画面座標へ写します。物理と表示を混ぜると、パラメータを変えたときに壊れます

const x1 = p.L1 * Math.sin(s.theta1)
const y1 = p.L1 * Math.cos(s.theta1)
const x2 = x1 + p.L2 * Math.sin(s.theta2)
const y2 = y1 + p.L2 * Math.cos(s.theta2)

鉛直下向きを θ = 0 に取ると、y = L cos θ で下向きが正になります。画面の y 軸も下向きが正なので、この取り方だと符号の反転が要りません。地味ですが、こういう座標系の選択でバグが減ります。

#😓 苦労したところ

カオスであることをどう見せるか。 「初期値をわずかに変えると全然違う動きになる」が売りなのに、1 本だけ動かしても伝わりません。プリセットで「ほぼ同じ初期条件」を用意して、軌跡の違いで体感してもらう形にしました。

減衰を入れるか。 現実の振り子は摩擦で止まります。入れないと永久に動き続けて、軌跡が画面を埋め尽くす。減衰係数をパラメータにして、0 にもできるようにしました。ただし減衰を入れるとエネルギー保存は成り立たなくなるので、検算表示は減衰 0 のときだけ意味がある、と注記しています。

#🔭 今後の拡張

• 2 本の振り子を並べて、初期値のわずかな差が開いていく様子を見せる
• リアプノフ指数の推定(カオスの度合いを数値で出す)
• 三重振り子
• 軌跡の書き出し(SVG / 動画)

#💡 このサービスから言えること

数値計算を出すサービスは、検算の手段とセットで作るべきです。 二重振り子は「それらしく動いている」だけなら簡単に作れますが、それが正しいかは見た目では分かりません。エネルギー保存という独立した確認方法があったから、式の写し間違いに気づけました。

正しさを確かめる方法が用意できない題材なら、そもそも数字を出すサービスにしない ── というのが、ここから引いている線です。