shigureni.dev

重力で引き合う — N 体

質点が二つあって、互いに引き合う。引力の強さは距離の二乗に反比例する。F = G·m1·m2/r²。近いほど強く、遠いほど弱い。この力を加速度として速度に足し、速度を位置に足すと、二点が引かれて軌道を描く。

万有引力。質量を持つ点どうしに働く引力で、大きさは G·m1·m2/r²。G は引力の強さを決める定数、r は二点間の距離。各点に働く力は、自分以外の全点からの引力のベクトル和。点が三つ以上だと解が閉じた式で書けない（三体問題）。力を足して数値で積分する。

二点を左右に置いて、互いに垂直な初速を与える。毎フレーム相手への引力を求めて、互いに逆向きの加速度として速度に足す。位置に速度を足すと、二点が引かれながら回る。

export default (c) => {
  const w = c.canvas.width
  const h = c.canvas.height
  const cx = w / 2
  const cy = h / 2
  const G = 60
  const soft = 24
  const dt = 0.6

const bodies = [
    { x: cx - 90, y: cy, vx: 0, vy: -1.1, m: 3 },
    { x: cx + 90, y: cy, vx: 0, vy: 1.1, m: 3 },
  ]

c.fillStyle = colors.paper
  c.fillRect(0, 0, w, h)

return () => {
    c.fillStyle = colors.tintAlpha(238, 0.06)
    c.fillRect(0, 0, w, h)

const a = bodies[0]
    const b = bodies[1]
    const dx = b.x - a.x
    const dy = b.y - a.y
    const d2 = dx * dx + dy * dy + soft * soft
    const inv = 1 / (d2 * Math.sqrt(d2))
    const fa = G * b.m * inv
    const fb = G * a.m * inv
    a.vx += dx * fa * dt
    a.vy += dy * fa * dt
    b.vx -= dx * fb * dt
    b.vy -= dy * fb * dt
    a.x += a.vx * dt
    a.y += a.vy * dt
    b.x += b.vx * dt
    b.y += b.vy * dt

for (const p of bodies) {
      c.fillStyle = colors.tint(28)
      c.beginPath()
      c.arc(p.x, p.y, 2 + p.m, 0, Math.PI * 2)
      c.fill()
    }
  }
}

dx, dy が相手への向き。d2 が距離の二乗で、1 / (d2 * Math.sqrt(d2)) が 1/r³。向き dx, dy に 1/r³ を掛けると、長さが 1/r² で向き付きの加速度になる（dx の長さが r なので、r/r³ = 1/r²）。a には相手 b の質量ぶん、b には a の質量ぶんを逆向きに足すので、二点は互いに引き合う。半透明の paper を毎フレーム重ねていて、通り道が尾を引く。vy の初速を上げると軌道がほどけて楕円が伸び、下げると正面衝突に近づく。

d2 に足している soft * soft が効くのは二点が接近したとき。

1/r² は r が 0 に近づくと無限に発散する。二点がほぼ重なった瞬間、力が爆発的に大きくなり、一フレームで速度が跳ね上がって軌道が壊れる（数値積分は連続でなく離散の刻みなので、すり抜けと暴発が起きる）。ソフトニング。r² に下限の定数 ε² を足した 1/(r²+ε²) で評価する手で、近距離で力が頭打ちになり、衝突しても発散しない。ε（softening length）より近い距離は潰れて区別されなくなる。

ソフトニングの効きを見るために、二点を正面からぶつける。soft を時間で 6 から 36 のあいだで往復させると、同じ衝突コースでも近距離での跳ね方が変わる。速度に上限を掛けて、刻みが粗くても有界に保つ。

export default (c) => {
  const w = c.canvas.width
  const h = c.canvas.height
  const cx = w / 2
  const cy = h / 2
  const G = 80
  const dt = 0.5
  let t = 0

const reset = () => [
    { x: cx - 120, y: cy - 4, vx: 1.0, vy: 0, m: 3 },
    { x: cx + 120, y: cy + 4, vx: -1.0, vy: 0, m: 3 },
  ]
  let bodies = reset()

c.fillStyle = colors.paper
  c.fillRect(0, 0, w, h)

return () => {
    t += 0.01
    const soft = 6 + (1 - Math.cos(t)) * 15

c.fillStyle = colors.tintAlpha(238, 0.05)
    c.fillRect(0, 0, w, h)

for (const p of bodies) {
      const sp = Math.hypot(p.vx, p.vy)
      if (sp > 6) {
        p.vx = (p.vx / sp) * 6
        p.vy = (p.vy / sp) * 6
      }
      p.x += p.vx * dt
      p.y += p.vy * dt
    }

if (
      a.x < -20 ||
      a.x > w + 20 ||
      b.x < -20 ||
      b.x > w + 20 ||
      a.y < -20 ||
      a.y > h + 20 ||
      b.y < -20 ||
      b.y > h + 20
    ) {
      bodies = reset()
    }

for (const p of bodies) {
      c.fillStyle = colors.tint(28)
      c.beginPath()
      c.arc(p.x, p.y, 2 + p.m, 0, Math.PI * 2)
      c.fill()
    }
  }
}

soft が小さいときは二点がすれ違う瞬間に強く弾かれて鋭く折れ、大きいときは近距離の力が鈍って緩く曲がる。画面外へ振り飛ばされたら初期配置に戻して、衝突を繰り返す。Math.hypot(p.vx, p.vy) で速さを測り、6 を越えたら向きを保ったまま長さを 6 に切り詰めている。刻み dt が粗いと一歩で進みすぎて発散するので、上限で頭を抑える。

点を二つから多数に増やす。各点について、自分以外の全点からの引力をベクトルで足し合わせる。N 個なら一フレームに N² 回の引力計算。中心へ弱く引く項を足して、群れが画面から散らないようにつなぎ止める。

export default (c) => {
  const w = c.canvas.width
  const h = c.canvas.height
  const cx = w / 2
  const cy = h / 2
  const N = 40
  const G = 22
  const soft = 26
  const dt = 0.6
  const cap = 4

const bodies = Array.from({ length: N }, () => {
    const ang = Math.random() * Math.PI * 2
    const r = 30 + Math.random() * (h * 0.28)
    return {
      x: cx + Math.cos(ang) * r,
      y: cy + Math.sin(ang) * r,
      vx: -Math.sin(ang) * (0.4 + Math.random() * 0.6),
      vy: Math.cos(ang) * (0.4 + Math.random() * 0.6),
      m: 0.6 + Math.random() * 2.4,
    }
  })

const accel = (i) => {
    const bi = bodies[i]
    let ax = (cx - bi.x) * 0.0009
    let ay = (cy - bi.y) * 0.0009
    for (let j = 0; j < N; j++) {
      if (j === i) continue
      const bj = bodies[j]
      const dx = bj.x - bi.x
      const dy = bj.y - bi.y
      const d2 = dx * dx + dy * dy + soft * soft
      const inv = 1 / (d2 * Math.sqrt(d2))
      const f = G * bj.m * inv
      ax += dx * f
      ay += dy * f
    }
    return [ax, ay]
  }

c.fillStyle = colors.paper
  c.fillRect(0, 0, w, h)

return () => {
    c.fillStyle = colors.tintAlpha(238, 0.07)
    c.fillRect(0, 0, w, h)

for (let i = 0; i < N; i++) {
      const b = bodies[i]
      const a = accel(i)
      b.vx += a[0] * dt
      b.vy += a[1] * dt
      const sp = Math.hypot(b.vx, b.vy)
      if (sp > cap) {
        b.vx = (b.vx / sp) * cap
        b.vy = (b.vy / sp) * cap
      }
      b.x += b.vx * dt
      b.y += b.vy * dt

const g = Math.max(8, 58 - b.m * 14) | 0
      c.fillStyle = colors.tint(g)
      c.beginPath()
      c.arc(b.x, b.y, 1 + b.m, 0, Math.PI * 2)
      c.fill()
    }
  }
}

各点は円周上に撒いて、接線方向の初速を与えてある。accel が一点に働く加速度で、内側のループで自分以外の N-1 個の引力を足している。質量 m の重い点ほど暗く大きく塗り、周りの軽い点を強く引き寄せる。(cx - bi.x) * 0.0009 の中心への弱い引きが、全体を画面の中に留めるおもり。重い点のまわりに軽い点が巻きついて、渦を巻く軌跡が積もる。質量の分布を変えると、引きつける核の数が変わる。