Алгоритм N-тіл Барнса-Гата з нуля

1. Задача N-тіл і складність

У гравітаційній симуляції N-тіл кожна частинка притягує кожну іншу частинку. Наївний алгоритм обчислює N(N−1)/2 парних взаємодій за крок:

2. Структура квадродерева / октодерева

Барнс-Гат розбиває простір рекурсивно. У 2D це квадродерево (4 нащадки); у 3D — це октодерево (8 нащадків). Кожен вузол представляє квадратну (або кубічну) область і є або:

• листком, що містить рівно одну частинку, або
• внутрішнім вузлом, нащадки якого покривають 4 (або 8) підобластей.

// Вузол 2D-квадродерева
class QNode {
  constructor(cx, cy, size) {
    this.cx = cx;        // x центру клітинки
    this.cy = cy;        // y центру клітинки
    this.size = size;    // піврозмір клітинки
    this.mass = 0;       // загальна маса в цій клітинці
    this.comx = 0;       // x центру мас
    this.comy = 0;       // y центру мас
    this.particle = null; // не null, якщо листок
    this.children = null; // null, якщо листок; [nw, ne, sw, se], якщо внутрішній
  }
}

3. Вставка в дерево

function insert(node, p) {
  if (node.mass === 0) {
    // Порожній вузол — зберігаємо частинку як листок
    node.particle = p;
    node.mass = p.mass;
    node.comx = p.x;
    node.comy = p.y;
    return;
  }

  if (node.children === null) {
    // Був листком — розбиваємо й перевставляємо наявну частинку
    subdivide(node);
    insertToChild(node, node.particle);
    node.particle = null;
  }

  // Вставляємо нову частинку у відповідного нащадка
  insertToChild(node, p);

  // Оновлюємо накопичений центр мас
  node.mass += p.mass;
  node.comx = (node.comx * (node.mass - p.mass) + p.x * p.mass) / node.mass;
  node.comy = (node.comy * (node.mass - p.mass) + p.y * p.mass) / node.mass;
}

function subdivide(node) {
  const h = node.size / 2;
  node.children = [
    new QNode(node.cx - h, node.cy + h, h),  // ПнЗ
    new QNode(node.cx + h, node.cy + h, h),  // ПнСх
    new QNode(node.cx - h, node.cy - h, h),  // ПдЗ
    new QNode(node.cx + h, node.cy - h, h),  // ПдСх
  ];
}

function insertToChild(node, p) {
  const isEast = p.x >= node.cx;
  const isNorth = p.y >= node.cy;
  const idx = (isNorth ? 0 : 2) + (isEast ? 1 : 0);
  insert(node.children[idx], p);
}

4. Обчислення центру мас

Положення центру мас і загальна маса, збережені в кожному внутрішньому вузлі, оновлюються під час вставки (див. вище). Вони представляють сукупний гравітаційний вплив усіх частинок у цій клітинці:

Коли ми обходимо дерево під час обчислення сил, віддалена клітинка трактується як єдина частинка з масою M_cell, розташована в (x_com, y_com). Це мультипольне наближення монопольного порядку — мультипольні методи вищого порядку (FMM) розширюють це квадрупольними та октупольними членами для більшої точності.

5. Критерій тета

Наскільки далеко — це «достатньо далеко», щоб застосувати наближення? Це вирішує критерій тета (критерій кута розкриття):

• θ = 0: ніколи не наближати — еквівалентно повному перебору O(N²).
• θ = 0.5: добрий баланс. Похибка ~1% для більшості симуляцій галактичного масштабу.
• θ = 1.0: швидше, але менш точно.
• θ = ∞: наближати все — використовувати лише центр мас кореня (жахливо).

6. Обхід для обчислення сил

const G = 6.674e-11;
const SOFTENING = 0.1;    // запобігає F → ∞ при r → 0
const THETA     = 0.5;

function calcForce(node, p, out) {
  if (node.mass === 0) return;

  const dx = node.comx - p.x;
  const dy = node.comy - p.y;
  const distSq = dx*dx + dy*dy + SOFTENING*SOFTENING;
  const dist   = Math.sqrt(distSq);

  if (node.children === null || node.size / dist < THETA) {
    // Листок АБО клітинка достатньо далеко — використовуємо наближення
    if (node.particle === p) return;  // пропускаємо саму себе
    const f = G * p.mass * node.mass / distSq;
    out.fx += f * dx / dist;
    out.fy += f * dy / dist;
  } else {
    // Клітинка надто близько — рекурсивно заходимо в нащадків
    for (const child of node.children) {
      if (child && child.mass > 0) calcForce(child, p, out);
    }
  }
}

7. Повний цикл інтегрування

function step(particles, dt) {
  // 1. Обчислюємо обмежувальний прямокутник усіх частинок
  let minX=Infinity, maxX=-Infinity, minY=Infinity, maxY=-Infinity;
  for (const p of particles) {
    if (p.x < minX) minX=p.x; if (p.x > maxX) maxX=p.x;
    if (p.y < minY) minY=p.y; if (p.y > maxY) maxY=p.y;
  }
  const cx   = (minX + maxX) / 2;
  const cy   = (minY + maxY) / 2;
  const size = Math.max(maxX - minX, maxY - minY) / 2 * 1.01;

  // 2. Будуємо квадродерево
  const root = new QNode(cx, cy, size);
  for (const p of particles) insert(root, p);

  // 3. Для кожної частинки обчислюємо силу та інтегруємо (leapfrog)
  for (const p of particles) {
    const force = { fx: 0, fy: 0 };
    calcForce(root, p, force);

    // Інтегрування методом Velocity Verlet / Leapfrog
    p.vx += (force.fx / p.mass) * dt;
    p.vy += (force.fy / p.mass) * dt;
    p.x  += p.vx * dt;
    p.y  += p.vy * dt;
  }
}

8. Поради з оптимізації

• Пулінг обʼєктів: виділіть плоский масив із MAX_NODES заздалегідь створених обʼєктів QNode і скидайте їх щокроку замість виділення мільйонів обʼєктів за секунду.
• Дерево у плоскому масиві: зберігайте дерево як плоский типізований масив (SoA — структура масивів) з цілочисловими індексами нащадків. Дружнє до кешу; уникає переходів за вказівниками на великих N.
• Web Workers: винесіть обчислення сил в окремий потік Worker (postMessage / SharedArrayBuffer). Цикл рендерингу в головному потоці читає результати, не блокуючи інтерфейс.
• GPU через обчислення WebGL: обчислювальні шейдери WebGPU можуть виконувати обхід BH повністю на GPU. Для WebGL 2 можливе сплощене дерево, збережене в текстурі та зчитуване обчислювально-проксі вершинним шейдером, хоч це й складно.
• Адаптивна тета: збільшуйте θ для частинок, що далеко від усіх центрів скупчень (де точність менш важлива), і зменшуйте θ у щільних областях.
• Помʼякшення: мала довжина помʼякшення ε (додана до знаменника як √(r² + ε²)) запобігає сингулярності при r=0 та уникає нереалістично великих прискорень між близькими частинками. Типове значення: 1–5% від середньої відстані між частинками.

З цими оптимізаціями однопотокова реалізація на JavaScript може опрацьовувати ~50 000 частинок з інтерактивною частотою кадрів, а реалізація з Web Worker може дійти до ~200 000.