Drum Modes: Laplacian Eigenfunctions

tap for next mode

The Dirichlet Laplacian $- Δ$ on a bounded planar domain $Ω$ has a discrete spectrum $0 < λ_{1} < λ_{2} \leq λ_{3} \leq \dots$ with eigenfunctions $ϕ_{k}$ vanishing on $\partial Ω$ :

- Δ ϕ_{k} = λ_{k} ϕ_{k} in Ω, ϕ_{k} = 0 on \partial Ω.

Physically these are the pure tones of a drumhead clamped along its rim. On the unit disk the eigenproblem separates in polar coordinates, and the eigenfunctions are exactly

ϕ_{n, m} (r, θ) = J_{n} (λ_{n, m} r) cos (n θ),

where

J_{n}

is the Bessel function of the first kind of order

n

and

λ_{n, m}

is its

m

-th positive zero. The corresponding eigenvalue is

λ_{n, m}^{2}

and the time evolution from the wave equation

u_{tt} = c^{2} Δ u

cos (ω_{n, m} t)

with

ω_{n, m} \propto λ_{n, m}

, so higher modes vibrate faster — that is the physics of drum overtones. The animation cycles through the first 16 modes ordered by

λ

; red is positive deflection, blue negative, and the cream curves are the nodal lines

{ϕ_{n, m} = 0}

(Courant's bound: at most

k

nodal domains for the

k

-th eigenfunction). Mark Kac asked in 1966 *can one hear the shape of a drum?* — the answer turned out to be no in general (Gordon–Webb–Wolpert, 1992), but the spectrum still encodes a great deal:

λ_{1}

scales as

1/ area

(Faber–Krahn), Weyl's law recovers area and perimeter from the counting function. Inspired by Gabriel Peyré's spectral-geometry visualizations.

idle

204 lines · vanilla

view source

// Drum modes — vibrating modes of a circular membrane (disk).
//
// The Dirichlet Laplacian on the unit disk has closed-form eigenfunctions
//   phi_{n,m}(r, theta) = J_n(lambda_{n,m} * r) * cos(n * theta)
// where lambda_{n,m} is the m-th positive zero of the Bessel function J_n.
// Eigenvalue: lambda_{n,m}^2. Time evolution: cos(omega t) with
// omega proportional to lambda_{n,m} (wave equation: u_tt = c^2 Delta u).
//
// Inspired by Gabriel Peyré's spectral-geometry visualizations.

// First 16 modes ordered by lambda (the eigenvalue spectrum of the unit disk).
// Each entry: [n, m, lambda_{n,m}]   (n = angular index, m = radial index,
// lambda = m-th positive zero of J_n).
const MODES = [
  [0, 1,  2.4048255577],
  [1, 1,  3.8317059702],
  [2, 1,  5.1356223019],
  [0, 2,  5.5200781103],
  [3, 1,  6.3801618959],
  [1, 2,  7.0155866698],
  [4, 1,  7.5883424345],
  [2, 2,  8.4172441404],
  [0, 3,  8.6537279129],
  [5, 1,  8.7714838159],
  [3, 2,  9.7610231300],
  [6, 1,  9.9361095242],
  [1, 3, 10.1734681351],
  [4, 2, 11.0647094885],
  [7, 1, 11.0863700450],
  [2, 3, 11.6198411721],
];

const MODE_DURATION = 3.0; // seconds per mode (before auto-advance)

let W, H;
let imgBuf;        // ImageData for the disk render (square, side = SIZE)
let SIZE;          // pixel side of the rendered disk image (recomputed on resize)
let modeIdx;
let modeElapsed;   // seconds in the current mode
let prevClickConsumed;

// Bessel function J_n(x) via series for small x, asymptotic for large x.
function J0(x) {
  const ax = Math.abs(x);
  if (ax < 8.0) {
    const y = x * x;
    // Abramowitz & Stegun 9.4.1
    const p = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
      + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
    const q = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
      + y * (59272.64853 + y * (267.8532712 + y * 1.0))));
    return p / q;
  }
  const z = 8.0 / ax;
  const y = z * z;
  const p = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
  const q = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5
    + y * (0.7621095161e-6 + y * (-0.934935152e-7))));
  return Math.sqrt(0.636619772 / ax) * (Math.cos(ax - 0.785398164) * p
    - z * Math.sin(ax - 0.785398164) * q);
}

function J1(x) {
  const ax = Math.abs(x);
  if (ax < 8.0) {
    const y = x * x;
    const p = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1
      + y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
    const q = 144725228442.0 + y * (2300535178.0 + y * (18583304.74
      + y * (99447.43394 + y * (376.9991397 + y * 1.0))));
    return p / q;
  }
  const z = 8.0 / ax;
  const y = z * z;
  const p = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
    + y * (0.2457520174e-5 + y * (-0.240337019e-6))));
  const q = 0.04687499995 + y * (-0.2002690873e-3 + y * (0.8449199096e-5
    + y * (-0.88228987e-6 + y * 0.105787412e-6)));
  let ans = Math.sqrt(0.636619772 / ax) * (Math.cos(ax - 2.356194491) * p
    - z * Math.sin(ax - 2.356194491) * q);
  if (x < 0) ans = -ans;
  return ans;
}

// J_n via upward recurrence is unstable for n > x; use Miller's downward
// recurrence for n >= 2.
function Jn(n, x) {
  if (n === 0) return J0(x);
  if (n === 1) return J1(x);
  if (x === 0) return 0;
  const ax = Math.abs(x);
  if (ax > n) {
    // Upward recurrence is stable here.
    let bjm = J0(ax);
    let bj = J1(ax);
    const tox = 2.0 / ax;
    for (let j = 1; j < n; j++) {
      const bjp = j * tox * bj - bjm;
      bjm = bj;
      bj = bjp;
    }
    return (x < 0 && (n & 1)) ? -bj : bj;
  }
  // Miller's downward recurrence.
  const tox = 2.0 / ax;
  const m = 2 * Math.floor((n + Math.floor(Math.sqrt(40 * n))) / 2);
  let jsum = 0;
  let bjp = 0;
  let bj = 1;
  let ans = 0;
  for (let j = m; j > 0; j--) {
    const bjm = j * tox * bj - bjp;
    bjp = bj;
    bj = bjm;
    if (Math.abs(bj) > 1e10) {
      bj *= 1e-10;
      bjp *= 1e-10;
      ans *= 1e-10;
      jsum *= 1e-10;
    }
    if (jsum) jsum += bj;
    jsum = !jsum && (j & 1) ? bj : jsum + bj;
    if (j === n) ans = bjp;
  }
  jsum = 2.0 * jsum - bj;
  ans /= jsum;
  return (x < 0 && (n & 1)) ? -ans : ans;
}

function colormap(t, out, off) {
  // t in [-1, 1]: blue (neg) -> cream (zero) -> red (pos).
  const a = Math.max(-1, Math.min(1, t));
  const m = Math.abs(a);
  const c0 = [244, 238, 228];
  const c1 = a >= 0 ? [196, 50, 64] : [48, 88, 178];
  out[off]     = Math.round(c0[0] + (c1[0] - c0[0]) * m);
  out[off + 1] = Math.round(c0[1] + (c1[1] - c0[1]) * m);
  out[off + 2] = Math.round(c0[2] + (c1[2] - c0[2]) * m);
  out[off + 3] = 255;
}

// Rebuild the disk image for the given mode + amplitude.
function renderDisk(n, lambda, amp) {
  const data = imgBuf.data;
  const half = SIZE * 0.5;
  // Slight inset so the rim sits inside the image.
  const Rpix = half - 1;
  for (let j = 0; j < SIZE; j++) {
    const dy = (j + 0.5) - half;
    for (let i = 0; i < SIZE; i++) {
      const off = (j * SIZE + i) * 4;
      const dx = (i + 0.5) - half;
      const rr = Math.sqrt(dx * dx + dy * dy);
      if (rr > Rpix) {
        // Outside disk — transparent slot for the bg.
        data[off]     = 12;
        data[off + 1] = 13;
        data[off + 2] = 18;
        data[off + 3] = 255;
        continue;
      }
      const r = rr / Rpix;            // r in [0, 1]
      const theta = Math.atan2(dy, dx);
      const val = Jn(n, lambda * r) * Math.cos(n * theta) * amp;
      colormap(val, data, off);
    }
  }
}

function init({ canvas, ctx, width, height, input }) {
  W = width;
  H = height;
  SIZE = Math.max(64, Math.floor(Math.min(W, H) * 0.86));
  imgBuf = new ImageData(SIZE, SIZE);
  modeIdx = 0;
  modeElapsed = 0;
  prevClickConsumed = false;

  ctx.fillStyle = '#0c0d12';
  ctx.fillRect(0, 0, W, H);
}

function advanceMode() {
  modeIdx = (modeIdx + 1) % MODES.length;
  modeElapsed = 0;
}

function tick({ ctx, dt, time, width, height, input }) {
  if (width !== W || height !== H) {
    W = width;
    H = height;
    const newSize = Math.max(64, Math.floor(Math.min(W, H) * 0.86));
    if (newSize !== SIZE) {
      SIZE = newSize;
      imgBuf = new ImageData(SIZE, SIZE);
    }
  }

  // Click handling — advance to next mode.
  if (input && typeof input.consumeClicks === 'function') {
    const clicks = input.consumeClicks();
    if (clicks && clicks > 0) advanceMode();
  } else if (input && input.mouseDown && !prevClickConsumed) {
    advanceMode();
    prevClickConsumed = true;
  } else if (input && !input.mouseDown) {
    prevClickConsumed = false;
  }

  modeElapsed += dt;
  if (modeElapsed >= MODE_DURATION) advanceMode();

  const [n, m, lambda] = MODES[modeIdx];
  // Time modulation: omega proportional to lambda (wave equation).
  // Scaled so mode 1 oscillates at a comfortable ~0.6 Hz.
  const omega = lambda * 0.55;
  const amp = Math.cos(omega * time);

  // Bg
  ctx.fillStyle = '#0c0d12';
  ctx.fillRect(0, 0, W, H);

  renderDisk(n, lambda, amp);

  // Paint the SIZE x SIZE image centered.
  const oc = new OffscreenCanvas(SIZE, SIZE);
  const octx = oc.getContext('2d');
  octx.putImageData(imgBuf, 0, 0);
  const x0 = (W - SIZE) * 0.5;
  const y0 = (H - SIZE) * 0.5;
  ctx.drawImage(oc, x0, y0);

  // Thin rim
  ctx.save();
  ctx.strokeStyle = 'rgba(220, 220, 230, 0.7)';
  ctx.lineWidth = 1.0;
  ctx.beginPath();
  ctx.arc(x0 + SIZE * 0.5, y0 + SIZE * 0.5, SIZE * 0.5 - 1, 0, Math.PI * 2);
  ctx.stroke();
  ctx.restore();

  // HUD
  ctx.fillStyle = 'rgba(232, 228, 218, 0.95)';
  ctx.font = '13px ui-monospace, SFMono-Regular, Menlo, monospace';
  ctx.textBaseline = 'top';
  ctx.fillText(`mode ${modeIdx + 1} / ${MODES.length}`, 14, 12);
  ctx.fillText(`(n, m) = (${n}, ${m})`, 14, 30);
  ctx.fillText(`λ = ${lambda.toFixed(4)}`, 14, 48);
  ctx.fillStyle = 'rgba(232, 228, 218, 0.55)';
  ctx.font = '11px ui-sans-serif, system-ui, sans-serif';
  ctx.fillText('Dirichlet −Δ on the unit disk · Jₙ(λr) cos(nθ) cos(ωt)', 14, H - 22);
}

Drum Modes: Laplacian Eigenfunctions

Comments (0)