RK4 у JavaScript: 30-рядковий інтегратор

1. Чому Ейлера недостатньо

Найпростіший спосіб чисельно проінтегрувати dy/dt = f(t, y) — це явний (прямий) метод Ейлера: взяти поточний нахил і зробити крок уздовж нього.

y += f(t, y) * dt;
t += dt;

Це один рядок, але локальна похибка становить O(dt²) на крок — метод використовує похідну лише на початку інтервалу, тож будь-яка кривизна справжньої траєкторії протягом цього кроку повністю пропускається. Для хаотичної системи на кшталт рівнянь Лоренца ця похибка не просто трохи спотворює траєкторію — оскільки близькі траєкторії розходяться експоненціально, невелика похибка інтегрування швидко перетворюється на зовсім іншу орбіту. При більшому кроці метод Ейлера також може бути нестабільним: енергія чи амплітуда дрейфують, доки симуляція візуально не «вибухає».

2. Чотири стадії RK4

Метод Рунге-Кутта 4-го порядку (RK4) обчислює похідну чотири рази за крок — на початку, у двох оцінках середини та в кінці — а потім комбінує їх із вагами 1 : 2 : 2 : 1:

Інтуїтивно: k1 — це той самий нахил, що дав би звичайний Ейлер. k2 та k3 уточнюють цю оцінку, обчислюючи нахил у середині інтервалу — використовуючи оцінку попередньої стадії, щоб туди дістатися. k4 обчислює нахил на дальньому кінці, використовуючи оцінку середини для проекції вперед. Усереднення всіх чотирьох — з подвійною вагою для оцінок середини — скасовує головні члени похибки, даючи локальну похибку O(dt⁵) на крок (глобальну — O(dt⁴)) лише за 4-кратну вартість одного кроку Ейлера.

3. 30-рядковий інтегратор

Секрет перевикористовуваної реалізації RK4 — тримати стан як звичайний масив і написати два крихітних векторних помічники (add, scale) замість того, щоб жорстко прописувати x, y, z. Тоді степер працює без змін, незалежно від того, чи y має 2 компоненти (маятник), 3 (Лоренц) чи 40 (ланцюжок зв'язаних осциляторів).

// --- векторні помічники -------------------------------------------------
function add(a, b) {
  return a.map((v, i) => v + b[i]);
}
function scale(a, s) {
  return a.map((v) => v * s);
}

// --- один крок RK4 -----------------------------------------------------
// f(t, y) -> масив похідних, тієї ж форми, що й y
function rk4Step(f, t, y, dt) {
  const k1 = f(t,              y);
  const k2 = f(t + dt / 2, add(y, scale(k1, dt / 2)));
  const k3 = f(t + dt / 2, add(y, scale(k2, dt / 2)));
  const k4 = f(t + dt,     add(y, scale(k3, dt)));

  const sum = k1.map(
    (_, i) => k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]
  );
  return add(y, scale(sum, dt / 6));
}

// --- запуск симуляції ---------------------------------------------
function integrate(f, y0, dt, steps) {
  let t = 0, y = y0, trail = [y0];
  for (let i = 0; i < steps; i++) {
    y = rk4Step(f, t, y, dt);
    t += dt;
    trail.push(y);
  }
  return trail;
}

Підрахуємо: два 3-рядкові векторні помічники, 11-рядковий rk4Step і 8-рядковий цикл-драйвер — рівно близько 30 рядків для повністю універсального ЗДР-інтегратора, без жодних залежностей.

4. Перевірка на системі Лоренца

Щоб проінтегрувати систему Лоренца, потрібна лише функція похідних, що відповідає сигнатурі f(t, y) вище — сам rk4Step не змінюється:

function lorenz(t, [x, y, z]) {
  const sigma = 10, rho = 28, beta = 8 / 3;
  return [
    sigma * (y - x),
    x * (rho - z) - y,
    x * y - beta * z,
  ];
}

const trail = integrate(lorenz, [0.1, 0, 0], 0.01, 20000);
// trail[i] = [x, y, z] на кроці i — подавайте прямо у WebGL line strip

Замініть lorenz на функцію похідних маятника ([θ, ω] → [ω, -g/L·sin θ]) чи систему Ресслера — та сама пара rk4Step/integrate продовжить працювати без змін — саме в цьому й полягає сенс зберігання стану як універсального масиву.

5. Вибір кроку інтегрування

Локальна похибка RK4 масштабується як dt⁵, тому зменшення dt удвічі скорочує похибку кроку приблизно до 1/32 попереднього значення — але подвоює кількість кроків, потрібних для покриття того самого проміжку часу. На практиці:

• dt = 0.005–0.01 добре працює для системи Лоренца при класичних параметрах (σ=10, ρ=28, β=8/3).
• Перевіряйте збереження, а не лише вигляд: для консервативних систем (маятник без загасання) побудуйте графік повної енергії у часі — якщо вона помітно дрейфує, зменшіть dt.
• Фіксований vs адаптивний крок: у цьому туторіалі використовується фіксований крок. Для промислової точності з меншою кількістю зайвих кроків адаптивний метод (наприклад, Дорман-Прінс, RK45) автоматично підлаштовує dt на основі оцінки локальної похибки — це надлишково для більшості інтерактивних WebGL-демо, але корисно знати про існування таких методів.

6. Типові пастки