Симуляція орбіти за 50 рядків JavaScript

1. Постановка задачі: нерухоме Сонце, вільна планета

Повна задача двох тіл дозволяє обом тілам рухатися навколо їхнього спільного центру мас. Але якщо одне тіло (Сонце) у тисячі разів важче за інше (планета чи космічний апарат), воно практично не рухається — тож ми можемо закріпити його в центрі координат і інтегрувати лише положення r та швидкість v легшого тіла. Це те саме спрощення, що використовується у статті про орбітальну механіку Кеплера, але тут ми інтегруємо чисельно замість аналітичного розв'язання рівняння Кеплера — підхід, що узагальнюється на орбіти, збурені іншими тілами.

const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');
canvas.width = canvas.height = 500;

const MU = 4000;               // гравітаційний параметр G*M, підібраний під пікселі
let pos = { x: 180, y: 0 };  // планета стартує у перигелії
let vel = { x: 0, y: 5.6 };  // початкова тангенціальна швидкість

2. Ньютонівська гравітація як прискорення

Закон тяжіння Ньютона дає прискорення, що завжди спрямоване від планети до Сонця, масштабоване обернено пропорційно до квадрата відстані:

У вигляді функції — вона приймає поточне положення й повертає вектор прискорення:

function accel(p) {
  const r2 = p.x * p.x + p.y * p.y;
  const r = Math.sqrt(r2);
  const f = -MU / (r2 * r);   // -μ / r³, застосовується до кожної компоненти нижче
  return { x: f * p.x, y: f * p.y };
}

3. Чому інтегрування Ейлера не працює

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

// НЕ РОБІТЬ так для орбіт:
const a = accel(pos);
vel.x += a.x * dt; vel.y += a.y * dt;
pos.x += vel.x * dt; pos.y += vel.y * dt;

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

4. Velocity-Verlet: стабільний інтегратор

Velocity-Verlet виправляє це, розділяючи оновлення швидкості на два напівкроки навколо оновлення позиції. Це все ще лише три додаткові рядки, але метод є симплектичним — він зберігає величину, дуже близьку до повної енергії, на довільній кількості кроків, тож замкнута орбіта залишається замкнутою замість того, щоб розкручуватись:

function step(dt) {
  const a0 = accel(pos);
  pos.x += vel.x * dt + 0.5 * a0.x * dt * dt;
  pos.y += vel.y * dt + 0.5 * a0.y * dt * dt;
  const a1 = accel(pos);           // повторне обчислення у НОВІЙ позиції
  vel.x += 0.5 * (a0.x + a1.x) * dt;
  vel.y += 0.5 * (a0.y + a1.y) * dt;
}

Ключова відмінність від Ейлера: прискорення обчислюється двічі за крок — раз у старій позиції (a0) і раз у новій (a1) — а оновлення швидкості усереднює обидва значення. Саме на цій зміні побудований галузевий стандарт — leapfrog-інтегратор N-тіл.

5. Малювання сліду

Щоб побачити, як еліпс сам себе накреслює, зберігайте останні кілька сотень позицій у масиві та малюйте їх з прозорістю, що згасає до найстарішої точки:

const trail = [];
function draw() {
  ctx.fillStyle = '#0a0e1a';
  ctx.fillRect(0, 0, canvas.width, canvas.height);
  ctx.save();
  ctx.translate(canvas.width / 2, canvas.height / 2);

  trail.push({ x: pos.x, y: pos.y });
  if (trail.length > 300) trail.shift();
  trail.forEach((t, i) => {
    ctx.fillStyle = `rgba(129,140,248,${i / trail.length})`;
    ctx.fillRect(t.x, t.y, 2, 2);
  });

  ctx.fillStyle = '#f5c842';               // Сонце
  ctx.beginPath(); ctx.arc(0, 0, 10, 0, 7); ctx.fill();
  ctx.fillStyle = '#4fa3e0';               // планета
  ctx.beginPath(); ctx.arc(pos.x, pos.y, 5, 0, 7); ctx.fill();
  ctx.restore();
}

6. Повний код (50 рядків)

Вставте це у HTML-файл з <canvas id="c"></canvas> — код працює самостійно, без залежностей:

const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');
canvas.width = canvas.height = 500;

const MU = 4000;
let pos = { x: 180, y: 0 };
let vel = { x: 0, y: 5.6 };
const trail = [];

function accel(p) {
  const r2 = p.x * p.x + p.y * p.y;
  const r = Math.sqrt(r2);
  const f = -MU / (r2 * r);
  return { x: f * p.x, y: f * p.y };
}

function step(dt) {
  const a0 = accel(pos);
  pos.x += vel.x * dt + 0.5 * a0.x * dt * dt;
  pos.y += vel.y * dt + 0.5 * a0.y * dt * dt;
  const a1 = accel(pos);
  vel.x += 0.5 * (a0.x + a1.x) * dt;
  vel.y += 0.5 * (a0.y + a1.y) * dt;
}

function draw() {
  ctx.fillStyle = '#0a0e1a';
  ctx.fillRect(0, 0, canvas.width, canvas.height);
  ctx.save();
  ctx.translate(canvas.width / 2, canvas.height / 2);

  trail.push({ x: pos.x, y: pos.y });
  if (trail.length > 300) trail.shift();
  trail.forEach((t, i) => {
    ctx.fillStyle = `rgba(129,140,248,${i / trail.length})`;
    ctx.fillRect(t.x, t.y, 2, 2);
  });
  ctx.fillStyle = '#f5c842';
  ctx.beginPath(); ctx.arc(0, 0, 10, 0, 7); ctx.fill();
  ctx.fillStyle = '#4fa3e0';
  ctx.beginPath(); ctx.arc(pos.x, pos.y, 5, 0, 7); ctx.fill();
  ctx.restore();
}

function loop() {
  step(0.06);
  draw();
  requestAnimationFrame(loop);
}
loop();