Van der Pol Oscillator

// Van der Pol oscillator: x'' - mu(1-x^2) x' + x = 0.
// Top pane: phase portrait (x vs x') with fading trajectory spiraling
// onto the limit cycle. Bottom pane: time-domain x(t).
// Cursor X scrubs mu from 0.1 (near-sinusoidal) to 5 (relaxation).
const SUBSTEPS = 8;
const TRAIL_LEN = 1400;
const TS_LEN = 600;
const MU_MIN = 0.1;
const MU_MAX = 5.0;

let x, v, mu, t, trailX, trailV, head, count, tsX, tsHead, tsCount;

function reset(initX, initV) {
  x = initX; v = initV; t = 0;
  head = 0; count = 0;
  tsHead = 0; tsCount = 0;
  for (let i = 0; i < TRAIL_LEN; i++) { trailX[i] = x; trailV[i] = v; }
  for (let i = 0; i < TS_LEN; i++) tsX[i] = x;
}

function init() {
  trailX = new Float32Array(TRAIL_LEN);
  trailV = new Float32Array(TRAIL_LEN);
  tsX = new Float32Array(TS_LEN);
  mu = 1.0;
  reset(2.0, 0.0);
}

function deriv(xx, vv, m) {
  return [vv, m * (1 - xx * xx) * vv - xx];
}

function rk4(dt) {
  const m = mu;
  const k1 = deriv(x, v, m);
  const k2 = deriv(x + 0.5 * dt * k1[0], v + 0.5 * dt * k1[1], m);
  const k3 = deriv(x + 0.5 * dt * k2[0], v + 0.5 * dt * k2[1], m);
  const k4 = deriv(x + dt * k3[0], v + dt * k3[1], m);
  x += (dt / 6) * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]);
  v += (dt / 6) * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]);
  t += dt;
}

function tick({ ctx, dt, width, height, input }) {
  // Cursor X -> mu (continuous). Defaults to mu=1 before first hover.
  if (input.mouseX >= 0 && input.mouseX <= width) {
    const f = Math.max(0, Math.min(1, input.mouseX / width));
    mu = MU_MIN + f * (MU_MAX - MU_MIN);
  }
  if (input.consumeClicks().length > 0) reset(2.0, 0.0);

  // Integrate. Smaller step for stiff (large mu) regime.
  const stepDt = Math.min(dt, 1 / 30);
  const sub = mu > 2.5 ? SUBSTEPS * 3 : SUBSTEPS;
  const h = stepDt / sub;
  for (let s = 0; s < sub; s++) {
    rk4(h);
    trailX[head] = x; trailV[head] = v;
    head = (head + 1) % TRAIL_LEN;
    if (count < TRAIL_LEN) count++;
  }
  tsX[tsHead] = x;
  tsHead = (tsHead + 1) % TS_LEN;
  if (tsCount < TS_LEN) tsCount++;

  // Background.
  ctx.fillStyle = "#05080f";
  ctx.fillRect(0, 0, width, height);

  // Layout: phase portrait on top (~70%), time series below.
  const phaseH = Math.floor(height * 0.7);
  const tsY0 = phaseH + 8;
  const tsH = height - tsY0 - 22;

  // Phase portrait axes. v amplitude grows with mu; pick a range that
  // contains the limit cycle for all mu in our range.
  const vRange = 1.2 + mu * 1.3;          // ~[1.3, 7.7]
  const xRange = 2.6;                      // limit cycle amp ~2 for all mu
  const padL = 36, padR = 12, padT = 12, padB = 8;
  const plotX0 = padL, plotY0 = padT;
  const plotW = width - padL - padR;
  const plotH = phaseH - padT - padB;
  const cx = plotX0 + plotW / 2;
  const cy = plotY0 + plotH / 2;
  const sx = (plotW / 2) / xRange;
  const sy = (plotH / 2) / vRange;

  // Grid.
  ctx.strokeStyle = "rgba(120,140,180,0.10)";
  ctx.lineWidth = 1;
  ctx.beginPath();
  for (let gx = -2; gx <= 2; gx++) {
    const px = cx + gx * sx;
    ctx.moveTo(px, plotY0); ctx.lineTo(px, plotY0 + plotH);
  }
  const vStep = vRange > 4 ? 2 : 1;
  for (let gv = -Math.floor(vRange); gv <= Math.floor(vRange); gv += vStep) {
    const py = cy - gv * sy;
    ctx.moveTo(plotX0, py); ctx.lineTo(plotX0 + plotW, py);
  }
  ctx.stroke();

  // Axes.
  ctx.strokeStyle = "rgba(200,215,240,0.35)";
  ctx.lineWidth = 1.2;
  ctx.beginPath();
  ctx.moveTo(plotX0, cy); ctx.lineTo(plotX0 + plotW, cy);
  ctx.moveTo(cx, plotY0); ctx.lineTo(cx, plotY0 + plotH);
  ctx.stroke();

  ctx.fillStyle = "rgba(180,195,225,0.7)";
  ctx.font = "11px monospace";
  ctx.fillText("x'", cx + 4, plotY0 + 12);
  ctx.fillText("x", plotX0 + plotW - 10, cy - 4);

  // Trajectory: oldest -> newest with fading alpha and hue ramp.
  // Drawn as connected segments so the curve stays smooth.
  if (count > 1) {
    const start = (head - count + TRAIL_LEN) % TRAIL_LEN;
    ctx.lineWidth = 1.3;
    let prevPx = 0, prevPy = 0, havePrev = false;
    for (let i = 0; i < count; i++) {
      const idx = (start + i) % TRAIL_LEN;
      const px = cx + trailX[idx] * sx;
      const py = cy - trailV[idx] * sy;
      if (havePrev) {
        const a = 0.05 + 0.85 * (i / count);
        const hue = 195 + (i / count) * 80;
        ctx.strokeStyle = `hsla(${hue}, 90%, 65%, ${a})`;
        ctx.beginPath();
        ctx.moveTo(prevPx, prevPy);
        ctx.lineTo(px, py);
        ctx.stroke();
      }
      prevPx = px; prevPy = py; havePrev = true;
    }
  }

  // Head dot.
  const hx = cx + x * sx, hy = cy - v * sy;
  const glow = ctx.createRadialGradient(hx, hy, 0, hx, hy, 14);
  glow.addColorStop(0, "rgba(255,240,200,0.95)");
  glow.addColorStop(1, "rgba(255,200,120,0)");
  ctx.fillStyle = glow;
  ctx.beginPath(); ctx.arc(hx, hy, 14, 0, Math.PI * 2); ctx.fill();
  ctx.fillStyle = "#fff7d6";
  ctx.beginPath(); ctx.arc(hx, hy, 3.2, 0, Math.PI * 2); ctx.fill();

  // Time-series strip.
  ctx.fillStyle = "rgba(255,255,255,0.025)";
  ctx.fillRect(padL, tsY0, plotW, tsH);
  ctx.strokeStyle = "rgba(255,255,255,0.10)";
  ctx.strokeRect(padL + 0.5, tsY0 + 0.5, plotW - 1, tsH - 1);
  ctx.strokeStyle = "rgba(255,255,255,0.12)";
  ctx.beginPath();
  ctx.moveTo(padL, tsY0 + tsH / 2); ctx.lineTo(padL + plotW, tsY0 + tsH / 2);
  ctx.stroke();

  if (tsCount > 1) {
    ctx.strokeStyle = "hsl(200, 90%, 70%)";
    ctx.lineWidth = 1.4;
    ctx.beginPath();
    const startTs = (tsHead - tsCount + TS_LEN) % TS_LEN;
    for (let i = 0; i < tsCount; i++) {
      const idx = (startTs + i) % TS_LEN;
      const xx = padL + (i / (TS_LEN - 1)) * plotW;
      const norm = Math.max(-1, Math.min(1, tsX[idx] / 2.5));
      const yy = tsY0 + tsH / 2 - norm * (tsH / 2 - 2);
      if (i === 0) ctx.moveTo(xx, yy); else ctx.lineTo(xx, yy);
    }
    ctx.stroke();
  }

  ctx.fillStyle = "rgba(200,215,240,0.65)";
  ctx.font = "10px monospace";
  ctx.fillText("x(t)", padL + 4, tsY0 + 11);

  // HUD.
  ctx.fillStyle = "rgba(225,235,250,0.85)";
  ctx.font = "12px monospace";
  const regime = mu < 0.4 ? "near-sinusoidal" : mu < 1.6 ? "weakly nonlinear" : mu < 3 ? "nonlinear" : "relaxation";
  ctx.fillText(
    `mu=${mu.toFixed(2)}  ${regime}  t=${t.toFixed(1)}s  (move cursor; click to reset)`,
    12, height - 6
  );
}