Двовимірне хвильове рівняння — скінченні різниці FTCS, ping-pong у WebGL і рендеринг із картами нормалей

1. Хвильове рівняння — виведення та інтуїція

Двовимірне скалярне хвильове рівняння керує полем висот u(x, y, t) — тим, як скалярна величина (тиск, зміщення, електричне поле) еволюціонує в просторі та часі:

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

Загальний тривимірний розв’язок (Д’Аламбер) — це будь-яка суперпозиція функцій вигляду f(x − ct) та g(x + ct) (хвилі, що рухаються ліворуч і праворуч). У двовимірному випадку ситуація багатша — хвилі можуть поширюватися в усіх напрямках і інтерферувати, відбиватися та дифрагувати.

2. Скінченно-різницева схема FTCS — явне правило оновлення

Дискретизуємо простір із кроком сітки Δx і час із кроком Δt. Нехай uⁿᵢⱼ = u(iΔx, jΔy, nΔt). Замінюючи кожну часткову похідну її центрально-різницевою апроксимацією другого порядку:

Це класичний шаблон FTCS (Forward Time, Centred Space — вперед у часі, центровано в просторі). Оновлення кожної комірки залежить лише від її двох попередніх часових кроків і чотирьох безпосередніх сусідів — 5-точковий шаблон, що чудово відображається на семплер текстур GLSL із чотирма зчитуваннями зі зсувом.

3. Стійкість — умова Куранта–Фрідріхса–Леві

Схема FTCS умовно стійка. Умова CFL (Куранта–Фрідріхса–Леві) вимагає, щоб хвиля проходила щонайбільше одну комірку сітки за часовий крок:

Порушення умови CFL призводить до експоненційного зростання за кілька часових кроків — симуляція «вибухає». На практиці додають невеликий запас: беруть r ≈ 0.5.

4. Граничні умови

Поведінка на краях сітки визначає, як хвилі відбиваються або залишають область:

• Дiрixлe (закріплена): u = 0 на всіх граничних комірках. Хвилі відбиваються з інверсією фази, як струна, прив’язана до стіни. Просто реалізувати — достатньо щокроку обнуляти крайові комірки.
• Неймана (вільна): ∂u/∂n = 0 на межі — нульова нормальна похідна. Відображаємо значення внутрішньої комірки назовні (дзеркальне доповнення). Хвилі відбиваються у фазі, як вільний кінець струни.
• Періодична: загортання — правий край є сусідом лівого краю. Корисно для симуляцій нескінченної площини та розв’язувачів на основі ШПФ.
• Поглинальна (буферний шар): множимо u на коефіцієнт загасання, що плавно зростає від 1 усередині до 0 на краю в межах перехідної смуги. Мінімізує штучні відбиття — необхідно для моделювання відкритої води.

// Поглинальний буферний шар — попередньо обчислюємо вагу загасання для кожної комірки
function buildSponge(w, h, margin) {
  const d = new Float32Array(w * h);
  for (let y = 0; y < h; y++) {
    for (let x = 0; x < w; x++) {
      const ex = Math.min(x, w-1-x) / margin;   // 0..1 поблизу краю
      const ey = Math.min(y, h-1-y) / margin;
      const t  = Math.min(Math.min(ex, ey), 1.0);
      d[y*w+x] = t * t * (3 - 2*t);  // плавний перехід smoothstep
    }
  }
  return d;
}

5. Реалізація на Canvas 2D

Найпростіший підхід використовує три сітки Float32Array (поточна, попередня, наступна) і буфер ImageData для відображення. Навіть за роздільної здатності 256×256 це дає плавні брижі на 60 кадр/с.

const W = 256, H = 256;
const R = 0.5;   // квадрат числа Куранта = (c·dt/dx)²
const DAMP = 0.995;

let cur  = new Float32Array(W * H);
let prev = new Float32Array(W * H);
let next = new Float32Array(W * H);

function step() {
  for (let y = 1; y < H-1; y++) {
    for (let x = 1; x < W-1; x++) {
      const i = y * W + x;
      const laplacian =
        cur[i-1] + cur[i+1] + cur[i-W] + cur[i+W] - 4 * cur[i];
      next[i] = (2 * cur[i] - prev[i] + R * R * laplacian) * DAMP;
    }
  }
  // Циклічно міняємо буфери без виділення пам’яті
  [prev, cur, next] = [cur, next, prev];
}

function render(ctx) {
  const img = ctx.createImageData(W, H);
  const d   = img.data;
  for (let i = 0, n = W*H; i < n; i++) {
    const v = (cur[i] * 127 + 128) | 0;  // відображення [-1,1] → [1,255]
    d[i*4]   = 30 + v * 0.6 | 0;
    d[i*4+1] = 90 + v * 0.4 | 0;
    d[i*4+2] = v;
    d[i*4+3] = 255;
  }
  ctx.putImageData(img, 0, 0);
}

// Кидаємо ґаусів імпульс у позицію (px, py)
function addDrop(px, py, amp = 1.0, sigma = 5) {
  for (let y = py-15; y <= py+15; y++) {
    for (let x = px-15; x <= px+15; x++) {
      if (x < 0 || y < 0 || x >= W || y >= H) continue;
      const d2 = (x-px)*(x-px) + (y-py)*(y-py);
      cur[y*W+x] += amp * Math.exp(-d2 / (2*sigma*sigma));
    }
  }
}

За 256×256 × 60 кадр/с це ≈ 4 мільйони операцій множення-додавання на секунду — легко вкладається в бюджет CPU. Масштабуйте до 512×512, і JS починає не встигати; ось тоді стає потрібним шлях через GPU.

6. Ping-pong-обчислення у WebGL2

Для великих роздільних здатностей (512×512 і вище) ми переносимо оновлення хвилі на GPU за допомогою ping-pong-рендерингу в текстуру. Дві текстури з плаваючою комою зберігають поточне й попереднє поля висот; ми рендеримо в третій (або обмінюваний) фреймбуфер за допомогою фрагментного шейдера, що реалізує шаблон FTCS.

// Налаштування ping-pong у WebGL2
const gl = canvas.getContext('webgl2');

function makeRT(w, h) {
  const tex = gl.createTexture();
  gl.bindTexture(gl.TEXTURE_2D, tex);
  gl.texImage2D(gl.TEXTURE_2D, 0, gl.R32F, w, h, 0,
                 gl.RED, gl.FLOAT, null);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
  const fbo = gl.createFramebuffer();
  gl.bindFramebuffer(gl.FRAMEBUFFER, fbo);
  gl.framebufferTexture2D(gl.FRAMEBUFFER,
    gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, tex, 0);
  return { tex, fbo };
}

const rtA = makeRT(W, H);   // «поточний» буфер
const rtB = makeRT(W, H);   // «попередній» буфер
const rtC = makeRT(W, H);   // «наступний» вихідний буфер

Фрагментний шейдер оновлення — це шаблон FTCS, перекладений на GLSL:

// Фрагментний шейдер оновлення хвилі на GLSL (WebGL2, GLSL 300 es)
#version 300 es
precision highp float;
uniform sampler2D uCur, uPrev;
uniform float uR2;       // r² = (c·dt/dx)²
uniform float uDamp;
uniform vec2  uTexelSize; // vec2(1.0/W, 1.0/H) — розмір текселя
in vec2 vUv;
layout(location=0) out float fragHeight;

void main() {
  vec2 d = uTexelSize;
  float cur  = texture(uCur,  vUv).r;
  float prev = texture(uPrev, vUv).r;
  float L    = texture(uCur, vUv + vec2( d.x, 0.0)).r
             + texture(uCur, vUv + vec2(-d.x, 0.0)).r
             + texture(uCur, vUv + vec2(0.0,  d.y)).r
             + texture(uCur, vUv + vec2(0.0, -d.y)).r
             - 4.0 * cur;
  fragHeight = (2.0 * cur - prev + uR2 * L) * uDamp;
}

На кожному кадрі рендеримо повноекранний чотирикутник цим шейдером із текстурами A та B у C. Потім обертаємо: A ← B, B ← C. GPU виконує всі оновлення пікселів паралельно, масштабуючись до 2048×2048 на 60 кадр/с без участі CPU.

7. Генерація карт нормалей для 3D-рендерингу

Щоб відрендерити поле висот як освітлену 3D-поверхню (замість плоскої кольорової карти), обчисліть нормаль до поверхні в кожній точці. Для поля висот h(x, y) геометрична нормаль (до нормування) дорівнює:

Карту нормалей можна генерувати в другому проході фрагментного шейдера на кожному кадрі, зчитуючи поточну текстуру висот. У поєднанні з моделлю освітлення Фонґа чи Блінна–Фонґа й картою оточення для відбиттів ви отримуєте переконливу воду в реальному часі.

// Фрагментний шейдер нормалей + освітлення
#version 300 es
precision highp float;
uniform sampler2D uHeight;
uniform vec2  uTexelSize;
uniform vec3  uLightDir;
in vec2 vUv;
out vec4 fragColor;

void main() {
  vec2 d = uTexelSize;
  float hL = texture(uHeight, vUv - vec2(d.x, 0.0)).r;
  float hR = texture(uHeight, vUv + vec2(d.x, 0.0)).r;
  float hD = texture(uHeight, vUv - vec2(0.0, d.y)).r;
  float hU = texture(uHeight, vUv + vec2(0.0, d.y)).r;

  vec3 N = normalize(vec3(hL - hR, hD - hU, 0.04));
  float diff = max(dot(N, normalize(uLightDir)), 0.0);

  // Колір води: темно-синій + дзеркальний відблиск
  vec3 base = vec3(0.04, 0.18, 0.38);
  vec3 col  = base * (0.3 + 0.7*diff);
  vec3 H    = normalize(uLightDir + vec3(0.0, 0.0, 1.0));
  col += vec3(1.0) * pow(max(dot(N, H), 0.0), 48.0) * 0.8;

  fragColor = vec4(col, 1.0);
}

8. Розширення: загасання, дисперсія та зв’язок зі Шредінгером

Хвильове рівняння із загасанням

Додавання часової похідної першого порядку моделює в’язке загасання (втрату енергії хвилі):

На практиці для симуляцій у реальному часі достатньо щокроку множити все поле на скаляр, трохи менший за 1 (× 0.995), що моделює рівномірне в’язке затухання без зайвого ділення.

Дисперсійні хвилі

Реальні хвилі на воді дисперсійні — компоненти вищої частоти рухаються швидше за компоненти нижчої частоти. Дисперсійне співвідношення ω² = gk + σk³/ρ (гравітація + поверхневий натяг) не можна відтворити простим хвильовим рівнянням без додавання членів із вищими просторовими похідними. Поширене наближення — це рівняння Буссінеска або частотний підхід із застосуванням ШПФ (як у моделях океану Герстнера / IFFT).

Зв’язок із рівнянням Шредінгера

Замініть дійсне поле висот комплекснозначною хвильовою функцією ψ(x, y, t) і змініть ∂²/∂t² на i·∂/∂t:

TDSE використовує точно той самий просторовий лапласіанів шаблон, але еволюціонує комплексну фазу, а не дійсну амплітуду. Метод розщеплення оператора, що чергує напівкроки в координатному та імпульсному просторах (через ШПФ), дає симплектичний розв’язувач зі збереженням енергії — ідеальний для візуалізації квантового тунелювання, інтерференції на двох щілинах і дисперсії хвильового пакета.