Simulating the String

The wave equation tells you exactly how a plucked string moves. Turning that equation into code takes one key step: replace infinitely small differences with finite-sized ones. Everything else follows automatically.

The wave equation

A string under tension obeys Newton's second law at every point along it. The result is the 1D wave equation:

// ∂²u/∂t²  =  c²  ·  ∂²u/∂x²
//  acceleration    =  c²  ×  curvature

Here u(x, t) is the displacement of the string at position x and time t, and c is the wave speed — set by the string's tension and density. The equation says: wherever the string curves upward, it accelerates upward; wherever it curves downward, it accelerates down. That restoring force is what makes waves travel.

Exact solution: any initial shape splits into two mirror copies that travel left and right at speed c, bounce off the fixed endpoints, and pass through each other forever. The simulation does the same thing, one tiny step at a time.

Discretising: from equation to code

Replace the continuous string with N equally spaced nodes. Replace each derivative with the ratio of a small difference to a small step size. The second time-derivative becomes:

// ∂²u/∂t² ≈ ( u[t+1][i] − 2·u[t][i] + u[t−1][i] ) / dt²

The second space-derivative (curvature) becomes:

// ∂²u/∂x² ≈ ( u[t][i+1] − 2·u[t][i] + u[t][i−1] ) / dx²

Substitute both into c²·∂²x = ∂²t, then rearrange to isolate the only unknown — the future position u[t+1][i]:

u[t+1][i]  =  2·u[t][i]  −  u[t−1][i]
            +  r²·( u[t][i−1]  −  2·u[t][i]  +  u[t][i+1] )

// r = c·dt/dx   (Courant number — must stay ≤ 1 for stability)
// The future depends on: 1 past node + 3 present neighbours

The Courant number r is the simulation's stability dial. If r > 1, tiny errors amplify each step and the simulation explodes. At r = 1, information travels exactly one grid cell per time step — the ideal case matching the exact wave equation solution.

Step through the stencil

The four-point pattern — 1 node from the past row, 3 from the present — is called the stencil. Each future node needs exactly these four neighbours. Step through the computation one node at a time and watch the stencil advance across the grid.

t−1 past row

t present row

active stencil

t+1 computed

t+1 pending

Running at audio rate

Run the same update loop 44 100 times per second. Read one node as the audio output. The inner loop — sweeping the stencil across all N nodes — produces one audio sample per iteration.

for (t = 0; t < totalSamples; t++) {
  for (i = 1; i < N−1; i++) {
    future[i] = 2·present[i] − past[i]
               + r²·(present[i−1] − 2·present[i] + present[i+1]);
  }
  output[t] = present[outputNode];  // tap one node as audio
  [past, present, future] = [present, future, past];  // rotate buffers
}

String shape — all N nodes, updated each frame

Output node position matters. Tapping the string near an endpoint (bridge position) produces a bright, trebly sound with strong high harmonics. Tapping near the centre damps the odd harmonics, giving a rounder, flute-like tone — exactly how guitarists change timbre by picking position.

Karplus–Strong: the same thing, faster

The FD simulation keeps three arrays of N values — past, present, future. But look at the update: to compute future[i] you only need past[i] and three present neighbours. If you sweep left to right and write the result back into the same array as you go, the "past" values for the nodes you've already passed become the "future" values for the ones ahead. One circular buffer is enough.

That single ring of N samples — with a read head that laps around once per fundamental period — is the Karplus-Strong delay line. The averaging filter (½·(current + next)) approximates the FD curvature term: a gentle low-pass that removes energy faster from high frequencies, just like a real string's stiffness damping.

Memory: FD needs 3 × N floats; Karplus-Strong needs 1 × N. The time dimension is folded into the circular buffer. The output is perceptually identical — and it runs fast enough for real-time synthesis at any sample rate.

Ring buffer — each dot = one delay-line sample

Output waveform (last 2048 samples)