2D Wave Equation — FTCS finite differences, WebGL ping-pong and normal-map rendering

1. The wave equation — derivation and intuition

The 2D scalar wave equation governs a height field u(x, y, t) — how a scalar quantity (pressure, displacement, electric field) evolves in space and time:

Physically this says: acceleration equals curvature times a constant. A region that is a local maximum accelerates downward; a local minimum accelerates upward. This restoring force is exactly what drives oscillatory wave motion.

The general 3D solution (d'Alembert) is any superposition of functions of the form f(x − ct) and g(x + ct) (waves travelling left and right). In 2D the situation is richer — waves can spread in all directions and interfere, reflect, and diffract.

2. FTCS finite-difference scheme — the explicit update rule

We discretise space with grid spacing Δx and time with step Δt. Let uⁿᵢⱼ = u(iΔx, jΔy, nΔt). Replacing each partial derivative with its second-order centred-difference approximation:

This is the classic FTCS (Forward Time, Centred Space) stencil. The update for each cell depends only on its own previous two time-steps and its four immediate neighbours — a 5-point stencil that maps beautifully to a GLSL texture sampler with four offset reads.

3. Stability — the Courant-Friedrichs-Lewy condition

The FTCS scheme is conditionally stable. The CFL (Courant-Friedrichs-Lewy) condition requires the wave to travel at most one grid cell per time step:

Violating the CFL condition causes exponential blow-up in a few time steps — the simulation explodes. In practice add a small fudge factor: use r ≈ 0.5.

4. Boundary conditions

The behaviour at the edges of the grid determines how waves reflect or leave:

• Dirichlet (fixed): u = 0 at all boundary cells. Waves reflect with phase inversion, like a string tied to a wall. Simple to implement — just clamp border cells to zero each step.
• Neumann (free): ∂u/∂n = 0 at the boundary — zero normal derivative. Reflect the interior cell value outward (mirror padding). Waves reflect in phase, like a free string end.
• Periodic: wrap around — the right edge is the neighbour of the left edge. Useful for infinite-plane simulations and FFT-based solvers.
• Absorbing (sponge layer): multiply u by a damping factor that ramps from 1 in the interior to 0 at the edge over a transition band. Minimises artificial reflections — essential for simulating open water.

// Absorbing sponge layer — precompute per-cell damping weight
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 near edge
      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 ramp
    }
  }
  return d;
}

5. Canvas 2D implementation

The simplest approach uses three Float32Array grids (current, previous, next) and an ImageData buffer for display. Even at 256×256 resolution this gives smooth ripples at 60 fps.

const W = 256, H = 256;
const R = 0.5;   // Courant number squared = (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;
    }
  }
  // Rotate buffers without allocation
  [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;  // map [-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);
}

// Drop a Gaussian impulse at position (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));
    }
  }
}

At 256×256 × 60 fps this is ≈ 4 million multiply-adds per second — easily within CPU budget. Scale to 512×512 and JS starts to struggle; that's when the GPU path becomes necessary.

6. WebGL2 ping-pong compute

For large resolutions (512×512 or above) we move the wave update to the GPU using render-to-texture ping-pong. Two floating-point textures store the current and previous height fields; we render into a third (or swap) framebuffer using a fragment shader implementing the FTCS stencil.

// WebGL2 ping-pong setup
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);   // "current" buffer
const rtB = makeRT(W, H);   // "previous" buffer
const rtC = makeRT(W, H);   // "next" output buffer

The update fragment shader is the FTCS stencil translated to GLSL:

// Wave update GLSL fragment shader (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;
}

Each frame, render the screen quad using this shader with textures from A and B into C. Then rotate: A ← B, B ← C. The GPU executes all pixel updates in parallel, scaling to 2048×2048 at 60 fps with no CPU involvement.

7. Normal-map generation for 3D rendering

To render the height field as a lit 3D surface (instead of a flat colour map), compute the surface normal at each point. For a height field h(x, y), the geometric normal (before normalisation) is:

The normal map can be generated in a second fragment shader pass each frame, reading from the current height texture. Combined with a Phong or BlinnPhong lighting model and an environment map for reflections, you get convincing real-time water.

// Normal + lighting fragment shader
#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);

  // Water colour: dark blue + specular
  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. Extensions: damping, dispersion, and the Schrödinger connection

Damped wave equation

Adding a first-order time derivative models viscous damping (wave energy loss):

In practice for real-time simulations it is enough to multiply the entire field by a scalar slightly below 1 each step (× 0.995), which models uniform viscous attenuation without the extra division.

Dispersive waves

Real water waves are dispersive — higher-frequency components travel faster than lower-frequency ones. The dispersion relation ω² = gk + σk³/ρ (gravity + surface tension) cannot be captured by the simple wave equation without adding higher spatial derivative terms. A common approximation is the Boussinesq equations or a frequency-domain approach using the FFT (as in Gerstner / IFFT ocean models).

Link to the Schrödinger equation

Replace the real height field with a complex-valued wave function ψ(x, y, t) and change ∂²/∂t² to i·∂/∂t:

The TDSE uses exactly the same spatial Laplacian stencil but evolves the complex phase rather than a real amplitude. A split-operator method alternating half-steps in position and momentum space (via FFT) gives a symplectic, energy-conserving solver — perfect for visualising quantum tunnelling, double-slit interference, and wave-packet dispersion.