二重振り子のカオスを数値で見る — RK4積分と初期値鋭敏性の実装

#なぜ二重振り子なのか

単振り子は高校物理で扱う。周期は長さと重力だけで決まり、初期条件を少し変えても振る舞いはほぼ同じ。予測可能で、退屈なほど従順だ。

ところが振り子をもう一本つなげた瞬間、世界が変わる。二重振り子は決定論的カオスの代表例であり、運動方程式は完全に決定論的（ランダム要素ゼロ）なのに、初期角度を0.001ラジアンずらしただけで数秒後にはまったく異なる軌跡を描く。

DoublePendulum コンポーネント（362行）を実装する過程で、カオス力学の核心に何度もぶつかった。この記事では、ラグランジュ力学から始めて、RK4積分、NaN対策、初期値鋭敏性の可視化まで、実装者の視点で整理する。

#ラグランジュ力学から運動方程式へ

#なぜニュートン力学ではなくラグランジュか

ニュートン力学（F = ma）は自由な質点には強いが、拘束のある系では厄介だ。振り子は「棒の長さが一定」という拘束条件を持つ。ニュートン的に解くと張力という未知数が増え、方程式が複雑になる。

ラグランジュ力学は「エネルギーから力を導く」アプローチだ。運動エネルギー T と位置エネルギー V の差 L = T - V（ラグランジアン）を定義し、オイラー＝ラグランジュ方程式に代入すれば、拘束条件は自動的に処理される。

#状態ベクトル

二重振り子の状態は4つの変数で完全に記述できる:

• θ₁ (a1): 上の振り子の角度
• θ₂ (a2): 下の振り子の角度
• ω₁ (w1): 上の振り子の角速度
• ω₂ (w2): 下の振り子の角速度

interface State {
  a1: number  // θ₁
  a2: number  // θ₂
  w1: number  // ω₁ = dθ₁/dt
  w2: number  // ω₂ = dθ₂/dt
}

位置 (x, y) は角度と棒の長さから一意に決まるので、状態には含めない。これが一般化座標の利点だ。

#角加速度の導出

ラグランジアンから導出した角加速度の式を、実際のコードで使った形で示す。Δa = θ₁ - θ₂、m = M₁ + M₂ とする:

dω₁/dt:

dω₁/dt = [-g·m·sin(θ₁) - M₂·g·sin(θ₁ - 2θ₂)·0.5
          - sin(Δa)·M₂·(ω₂²·L₂ + ω₁²·L₁·cos(Δa))]
         / [L₁·(m - M₂·cos²(Δa))]

dω₂/dt:

dω₂/dt = [sin(Δa)·(ω₁²·L₁·m + g·m·cos(θ₁)
          + ω₂²·L₂·M₂·cos(Δa))]
         / [(L₂/L₁)·L₁·(m - M₂·cos²(Δa))]

式は長いが、やっていることはニュートンの第2法則の回転版だ。分子は「系に作用するトルクの合計」、分母は「有効な慣性モーメント」に対応する。

実装ではこれを derivatives 関数にまとめている:

function derivatives(s: State, g: number): State {
  const { a1, a2, w1, w2 } = s
  const da = a1 - a2
  const sinDa = Math.sin(da)
  const cosDa = Math.cos(da)
  const m = M1 + M2

  const dw1 =
    (-g * m * Math.sin(a1) -
      M2 * g * Math.sin(a1 - 2 * a2) * 0.5 -
      sinDa * M2 * (w2 * w2 * L2_RATIO + w1 * w1 * L1_RATIO * cosDa)) /
    (L1_RATIO * (m - M2 * cosDa * cosDa))

  const dw2 =
    (sinDa * (w1 * w1 * L1_RATIO * m + g * m * Math.cos(a1)
      + w2 * w2 * L2_RATIO * M2 * cosDa)) /
    ((L2_RATIO / L1_RATIO) * L1_RATIO * (m - M2 * cosDa * cosDa))

  return { a1: w1, a2: w2, w1: dw1, w2: dw2 }
}

戻り値の a1: w1, a2: w2 に注目してほしい。角度の時間微分は角速度そのものだ。この関数は「状態の時間微分」を返しており、これが積分器への入力になる。

#RK4積分の実装

#なぜオイラー法では駄目か

最もシンプルな数値積分法はオイラー法だ:

x(t + dt) = x(t) + f(x, t) * dt

一見よさそうだが、オイラー法は1次精度しかない。時間刻み dt を半分にしても誤差は半分にしかならない。カオス系では小さな誤差が指数的に増幅されるため、オイラー法ではすぐにエネルギーが発散する。振り子が回転し続けたり、突然停止したりと、物理的にあり得ない挙動になる。

#Runge-Kutta 4次法（RK4）

RK4は「4回の関数評価で4次精度」を実現する。dt を半分にすると誤差は1/16になる。費用対効果に優れている。

アイデアは、時間ステップの始点・中点・終点で傾きを評価し、加重平均をとることだ:

1. k₁: 現在の状態での傾き
2. k₂: k₁で半ステップ進めた状態での傾き
3. k₃: k₂で半ステップ進めた状態での傾き
4. k₄: k₃でフルステップ進めた状態での傾き
5. 最終更新: (k₁ + 2k₂ + 2k₃ + k₄) / 6

function rk4Step(s: State, g: number, dt: number): State {
  const k1 = derivatives(s, g)
  const s2: State = {
    a1: s.a1 + k1.a1 * dt * 0.5,
    a2: s.a2 + k1.a2 * dt * 0.5,
    w1: s.w1 + k1.w1 * dt * 0.5,
    w2: s.w2 + k1.w2 * dt * 0.5,
  }
  const k2 = derivatives(s2, g)
  const s3: State = {
    a1: s.a1 + k2.a1 * dt * 0.5,
    a2: s.a2 + k2.a2 * dt * 0.5,
    w1: s.w1 + k2.w1 * dt * 0.5,
    w2: s.w2 + k2.w2 * dt * 0.5,
  }
  const k3 = derivatives(s3, g)
  const s4: State = {
    a1: s.a1 + k3.a1 * dt,
    a2: s.a2 + k3.a2 * dt,
    w1: s.w1 + k3.w1 * dt,
    w2: s.w2 + k3.w2 * dt,
  }
  const k4 = derivatives(s4, g)

  return {
    a1: s.a1 + (k1.a1 + 2*k2.a1 + 2*k3.a1 + k4.a1) * (dt / 6),
    a2: s.a2 + (k1.a2 + 2*k2.a2 + 2*k3.a2 + k4.a2) * (dt / 6),
    w1: s.w1 + (k1.w1 + 2*k2.w1 + 2*k3.w1 + k4.w1) * (dt / 6),
    w2: s.w2 + (k1.w2 + 2*k2.w2 + 2*k3.w2 + k4.w2) * (dt / 6),
  }
}

#サブステップによる安定性確保

requestAnimationFrame は約60fpsで呼ばれる。1フレームあたり約16.7msだが、物理シミュレーションにはこの時間刻みは大きすぎる。

実装では DT = 0.02、STEPS_PER_FRAME = 8 としている。1フレームあたり8回のRK4ステップを回すことで、実効的な時間刻みを十分小さく保つ:

const DT = 0.02
const STEPS_PER_FRAME = 8

// アニメーションループ内
for (let i = 0; i < STEPS_PER_FRAME; i++) {
  stateRef.current = rk4Step(stateRef.current, g, DT)
}

DT = 0.02 の選択理由は実験的なものだ。0.05 ではエネルギーのドリフトが目視で分かる。0.01 では安定するが、STEPS_PER_FRAME を増やす必要がありCPU負荷が上がる。0.02 × 8 = 0.16 の実効ステップが、60fpsのフレーム間隔（~0.0167秒）に対してちょうどよい「少し速い時間の流れ」を作り、視覚的にも気持ちのいい速度感になる。

#NaN との戦い

#問題

浮動小数点演算には有限の精度しかない。IEEE 754 の64ビット浮動小数点数は約15桁の有効数字を持つが、カオス系では誤差が指数的に蓄積する。

具体的には、角加速度の分母 m - M₂·cos²(Δa) が 0 に近づくとき、除算結果が巨大になる。これが数ステップ続くと Infinity が生まれ、Infinity を含む計算は NaN を生む。NaN はあらゆる演算に伝播し、Canvas の描画が完全に止まる。

#対策

実装ではシンプルな Number.isFinite() ガードを入れている:

const { a1, a2, w1, w2 } = stateRef.current

if (
  !Number.isFinite(a1) ||
  !Number.isFinite(a2) ||
  !Number.isFinite(w1) ||
  !Number.isFinite(w2)
) {
  const angles = randomAngles()
  stateRef.current = { ...angles, w1: 0, w2: 0 }
  trailRef.current = []
  return
}

NaN/Infinity を検出したら、ランダムな初期角度で状態をリセットし、軌跡バッファもクリアする。ユーザーから見ると「振り子が新しい初期条件で再スタートした」ように見える。

カオス系は数値的にも不安定だ。ガードなしでは、重力を大きくした状態で数分放置するだけで破綻する。この種の防御的プログラミングは、物理シミュレーションでは必須と言っていい。

Number.isFinite() は NaN と ±Infinity の両方を弾く。!Number.isNaN() だけでは Infinity を見逃すので不十分だ。

#初期値鋭敏性の可視化

カオスの定義的な特性は初期値鋭敏性（sensitive dependence on initial conditions）だ。リアプノフ指数が正であること、つまり近接した軌道が指数的に発散することを意味する。

#速度→色のマッピング

軌跡の色は第2質点の速度（角速度の大きさ）に対応させている:

function speedToHue(speed: number): number {
  const t = Math.min(1, speed / 15)
  return 240 - t * 240
}

• 低速（speed ≈ 0）→ hue = 240（青）
• 高速（speed ≥ 15）→ hue = 0（赤）

青から赤へのグラデーションは直感的に「冷たい→熱い」「遅い→速い」と対応し、軌跡を見るだけでエネルギーの分布が分かる。

#軌跡バッファ

軌跡は TrailPoint[] の循環バッファとして実装している:

interface TrailPoint {
  x: number
  y: number
  speed: number
}

const MAX_TRAIL = 2000

// 毎フレーム
trailRef.current.push({ x: x2, y: y2, speed })
while (trailRef.current.length > MAX_TRAIL) {
  trailRef.current.shift()
}

2000点は60fpsで約33秒分の履歴にあたる。これ以上長くするとメモリと描画コストが増えるが、短すぎるとカオスの美しいパターンが見えない。2000は実験的に見つけたバランス点だ。

描画時には古い点ほど透明にすることで、軌跡が自然にフェードアウトする:

for (let i = 1; i < trail.length; i++) {
  const age = i / trail.length
  const alpha = age * 0.7 + 0.02
  const hue = speedToHue(trail[i].speed)

  ctx.beginPath()
  ctx.moveTo(trail[i - 1].x, trail[i - 1].y)
  ctx.lineTo(trail[i].x, trail[i].y)
  ctx.strokeStyle = `hsla(\${hue}, 80%, 60%, \${alpha})`
  ctx.stroke()
}

age は 0（最古）→ 1（最新）の正規化値で、新しい点ほど不透明になる。最古の点でも alpha = 0.02 としてわずかに残すことで、軌跡が突然消えるのを防いでいる。

#「0.001ラジアンの差」

実際にやってみると、初期角度を0.001ラジアン（約0.057度）ずらしただけで、5〜10秒後にはまったく異なる軌跡になる。人間の目には同じ初期条件に見えるのに、結果は劇的に違う。これがカオスだ。

この実装では Reset ボタンを押すたびに Math.random() で初期角度を生成するので、2回続けて同じ軌跡が出ることは実質的にあり得ない。すべての軌跡はユニークであり、一期一会だ。

#パラメータと操作感

#定数の設計

const M1 = 10       // 上の振り子の質量
const M2 = 10       // 下の振り子の質量
const L1_RATIO = 0.28  // ビューポート比での棒の長さ
const L2_RATIO = 0.28
const MAX_TRAIL = 2000
const DT = 0.02
const STEPS_PER_FRAME = 8

等質量（M1 = M2 = 10） を選んだ理由は、カオス性が最も顕著になるからだ。上の質量が極端に大きいと下の振り子は単なる付属物になり、カオスが減る。等質量では両者が対等に影響し合い、最も複雑な動きが生まれる。

L1_RATIO = L2_RATIO = 0.28 はビューポートの Math.min(width, height) に対する比率だ。0.28 × 2 = 0.56 なので、二本の棒を伸ばしきっても画面の56%に収まる。振り子の先端がキャンバスからはみ出さない最大値として設定した。

#重力スライダー

重力は 0.5〜3.0 の範囲で調整できる:

• g = 0.5: ゆっくりとした優雅な動き。軌跡が大きな曲線を描く
• g = 1.5（デフォルト）: 地球の重力をスケーリングした標準値
• g = 3.0: 激しい動き。角速度が大きくなり、NaN のリスクも高まる

重力を変えると系のエネルギースケールが変わるため、カオスの「見え方」も変わる。低重力ではゆったりとした渦巻き、高重力では激しい鞭のようなスナップが見られる。

#ピボット位置

const pivotX = w / 2
const pivotY = h * 0.3

ピボットを上方30%に配置することで、振り子の可動域を最大化している。中央に置くと下半分しか使えず、上方に置きすぎると棒が短くなりすぎる。

#実装で学んだこと

#ラグランジュ力学の実用性

ラグランジュ力学は「エネルギーから力を導く」アプローチだ。拘束条件（棒の長さが一定）を自動的に処理してくれるため、ニュートン力学のように張力を個別に計算する必要がない。一般化座標（角度）で記述すれば、運動方程式の導出は機械的な作業になる。

逆に言えば、自由度の多い系（N体問題など）ではニュートン力学のほうがシンプルなこともある。道具は問題に合わせて選ぶべきだ。

#カオス系のデバッグ

カオス系には「正解の軌跡」がない。2つの実装が同じ初期条件で異なる結果を出しても、どちらが正しいのか判定できない（微小な丸め誤差の違いが増幅されるため）。

代わりに定性的な性質で検証する:

1. エネルギー保存: 摩擦なしなら全エネルギー（運動 + 位置）は一定のはず。RK4でもわずかにドリフトするが、オイラー法よりはるかにマシ
2. 対称性: 左右反転した初期条件は左右反転した軌跡を生むはず
3. 極限ケース: 小角度（θ ≈ 0）では単振り子に近い周期運動になるはず

実装ではエネルギーの数値的なモニタリングは省いたが、開発中は console.log で T + V の値を追跡し、ドリフトが許容範囲内であることを確認した。

#RK4の費用対効果

RK4は「4回の関数評価で4次精度」を実現する。シンプレクティック積分器（Verlet法など）はエネルギー保存に優れるが、二重振り子のような非分離型ハミルトニアンには直接適用しにくい。RK4は汎用的で実装が簡単、かつ十分な精度があり、ブラウザでの物理シミュレーションにはベストバランスだと感じる。

4次の「4」は、テイラー展開の4次の項まで正確に再現するという意味だ。つまり誤差は O(dt⁵) であり、dt = 0.02 なら 1ステップあたりの誤差は約 3.2 × 10⁻⁹ 程度。これを1フレーム8回、60fps で積み上げても、数分間はエネルギーが実用上保存される。

#Canvas API描画の工夫

軌跡の描画では、セグメントごとに beginPath() → stroke() を呼んでいる。一見非効率に見えるが、各セグメントの色（hue）と透明度（alpha）が異なるため、一つの Path にまとめることができない。

パフォーマンスが問題になる場合は、軌跡を WebGL のラインストリップに置き換えるか、オフスクリーン Canvas に描画してキャッシュする手法が考えられるが、2000セグメント程度なら Canvas 2D で十分高速だ。

#まとめ：二重振り子とカオスの実装

二重振り子は、たった362行で決定論的カオスを体験できる題材だ。ラグランジュ力学で運動方程式を導き、RK4で積分し、NaN ガードで数値的安定性を確保し、速度→色マッピングで初期値鋭敏性を可視化する。

カオスは「予測不能」という意味で恐ろしいが、その裏には完全に決定論的な方程式がある。方程式は美しく、軌跡も美しい。ただし、その美しさを2回同じように見ることはできない。