N-body gravity: from O(N²) to Barnes-Hut

Newton's law of gravity

Гравітаційна сила між двома матеріальними точками масами m₁ та m₂, розділеними відстанню r, описується формулою Ньютона:

У симуляціях часто використовують безрозмірні одиниці або нормалізують G = 1, щоб уникнути проблем із числовою точністю float32. Прискорення тіла i від решти N–1 тіл:

Навіть проста сумація вимагає O(N²) обчислень за крок — для 10 000 тіл це вже 10⁸ операцій, що виводить за межі реального часу.

Softening: avoiding the singularity

Якщо два тіла близько зближуються, знаменник r³ підходить до нуля, і прискорення стає нескінченним — числова катастрофа. Рішення: додати softening-параметр ε (зазвичай позначають ε або ε²):

Фізичний сенс: тіла поводяться як м'які об'єкти зі «скругленим» потенціалом. Типові значення ε — від 0.01 до 1.0 у безрозмірних одиницях, залежно від масштабу задачі.

Naïve O(N²) algorithm

Базовий алгоритм — подвійний цикл: для кожної пари (i, j) обчислити силу і додати до вектора прискорення обох тіл. Завдяки 3-му закону Ньютона (actio = reactio) достатньо перебирати лише унікальні пари j > i:

Leapfrog integrator

Для N-body задач симплектичний Leapfrog (Störmer-Verlet DKD) є стандартом: він зберігає енергію набагато краще за метод Ейлера і дешевший за RK4.

Key property: symplecticity — Leapfrog exactly preserves the symplectic structure of phase space. This means orbits neither unwind nor spiral inward even after millions of steps.

Oct-Tree and Barnes-Hut

Алгоритм Barnes-Hut (1986) — ключова ідея: далекі тіла можна апроксимувати одним центром мас, а не підсумовувати кожне окремо.

Building the Oct-Tree

3D простір рекурсивно ділиться на 8 октантів доти, доки кожна клітина містить не більше одного тіла. Одночасно в кожній внутрішній ноді зберігають:

• Total mass M of all bodies in the cell
• Position of the centre of mass r̄
• Cell size s (cube edge length)

Tree traversal and the θ criterion

Для кожного тіла i дерево обходять рекурсивно. Для поточної ноди з розміром s і відстанню d від тіла перевіряють умову:

Де θ ≈ 0.5–1.0 — параметр точності. При θ = 0 отримуємо точний O(N²) алгоритм. При θ = 1.0 — максимальна швидкість, але помітна похибка в щільних регіонах.

θ criterion and accuracy

Для браузерної симуляції в реальному часі θ = 0.7–1.0 дає непомітну для ока похибку й достатню продуктивність. Наукові гідродинамічні коди зазвичай використовують θ = 0.5–0.7.

Complexity and benchmarks

Для 100 000 тіл Barnes-Hut дозволяє утримувати 50+ FPS, тоді як O(N²) займає 5 секунд на крок — повний гальм. Побудова дерева додає ~O(N log N) накладних витрат, але вони окупаються вже від ~500 тіл.

JavaScript implementation

Naïve O(N²) with Float32Array

// Float32Array for CPU cache efficiency:
const pos = new Float32Array(N * 3);   // x, y, z
const vel = new Float32Array(N * 3);
const acc = new Float32Array(N * 3);
const mass = new Float32Array(N);

function computeForces() {
  acc.fill(0);
  const eps2 = EPS * EPS;
  for (let i = 0; i < N; i++) {
    const ix = i * 3, iy = ix + 1, iz = ix + 2;
    for (let j = i + 1; j < N; j++) {
      const jx = j * 3, jy = jx + 1, jz = jx + 2;
      const dx = pos[jx] - pos[ix];
      const dy = pos[jy] - pos[iy];
      const dz = pos[jz] - pos[iz];
      const r2 = dx*dx + dy*dy + dz*dz + eps2;
      const inv = G * mass[j] / (r2 * Math.sqrt(r2));  // 1/r³
      acc[ix] += dx * inv;  acc[iy] += dy * inv;  acc[iz] += dz * inv;
      // Newton 3: force j on i is opposite
      const invI = G * mass[i] / (r2 * Math.sqrt(r2));
      acc[jx] -= dx * invI; acc[jy] -= dy * invI; acc[jz] -= dz * invI;
    }
  }
}

function leapfrogStep(dt) {
  const hdt = dt * 0.5;
  // Kick: v += a * dt/2
  for (let k = 0; k < N * 3; k++) vel[k] += acc[k] * hdt;
  // Drift: x += v * dt
  for (let k = 0; k < N * 3; k++) pos[k] += vel[k] * dt;
  // New forces at new positions
  computeForces();
  // Kick: v += a * dt/2
  for (let k = 0; k < N * 3; k++) vel[k] += acc[k] * hdt;
}

Oct-Tree node

class OctNode {
  constructor(cx, cy, cz, halfSize) {
    this.cx = cx; this.cy = cy; this.cz = cz;  // cube centre
    this.halfSize = halfSize;
    this.mass = 0;
    this.comX = 0; this.comY = 0; this.comZ = 0;   // centre of mass
    this.bodyIdx = -1;    // -1 = internal node
    this.children = null; // null or Array(8)
  }

  insert(i) {
    if (this.mass === 0) {          // empty node
      this.bodyIdx = i;
      this._updateCOM(i);
      return;
    }
    if (this.children === null) {  // leaf node with 1 body
      this._subdivide();
      this._insertIn(this.bodyIdx); // relocate existing body
      this.bodyIdx = -1;
    }
    this._insertIn(i);
    this._updateCOM(i);
  }

  calcAcc(i, ax, ay, az, theta) {
    if (this.mass === 0 || this.bodyIdx === i) return;
    const dx = this.comX - pos[i*3];
    const dy = this.comY - pos[i*3+1];
    const dz = this.comZ - pos[i*3+2];
    const d2 = dx*dx + dy*dy + dz*dz + EPS*EPS;
    // Barnes-Hut criterion: s/d < θ
    if ((this.bodyIdx >= 0) || (4 * this.halfSize * this.halfSize) / d2 < theta * theta) {
      const inv = G * this.mass / (d2 * Math.sqrt(d2));
      ax[i] += dx * inv;
      ay[i] += dy * inv;
      az[i] += dz * inv;
    } else {
      for (const child of this.children) {
        if (child) child.calcAcc(i, ax, ay, az, theta);
      }
    }
  }
}

Extensions and GPU

WebGPU and compute shaders

Наступний крок після Barnes-Hut — перенесення обчислень на GPU через WebGPU Compute Shaders. Пряме O(N²) на GPU (тисячі ядер) обганяє CPU Barnes-Hut до N ≈ 50 000. Вище — більш складні GPU-BVH або FMM алгоритми.

Fast Multipole Method (FMM)

FMM — наступник Barnes-Hut із справжньою O(N) складністю. Замість однієї апроксимації (маса + COM), кожна нода зберігає мультипольний розклад, що дозволяє точніше і швидше оцінювати далекі взаємодії. Однак реалізація значно складніша.

Adaptive time stepping

Тіла в щільних скупченнях потребують малого dt, тоді як далекі — великого. Block-based адаптивний крок часу, популярний у Gadget-2 та аналогах, може ще у 10–50 разів пришвидшити симуляцію без втрати точності в ключових регіонах.