The Kinked String Was Fine, the Cut String Overshot by 9% — How the 1747 Wave-Equation Fight Ended Numerically

How you represent a solution decides what kind of error you get. Travelling waves carry the shape; a sum of modes rings at every corner.

In 1747, three men split over a single string#

d'Alembert wrote down the equation of a vibrating string in 1747. It was the first partial differential equation anyone had written. The fight that followed was not about the equation. It was about what its solutions were allowed to look like. Euler and Daniel Bernoulli each produced a different answer, and none of the three backed down for thirty years.

All three were right. The answers only diverge once you cut the series off at a finite number of terms. This post puts both representations into the same code and watches how the smoothness of the initial data decides the kind of error you get. The 9% ringing you meet in spectral and high-order schemes starts here.

d'Alembert expanded nothing#

Apply Newton's second law to a small segment of a string with tension $T$ and line density $\rho$ and you get this.

u_{tt} = c^2 u_{xx}, \qquad c = \sqrt{T/\rho}

$u$ is the transverse displacement, $c$ the wave speed. d'Alembert noticed that the substitution $\xi = x - ct$ , $\eta = x + ct$ turns the equation into $u_{\xi\eta} = 0$ . Integrate twice and the solution falls out.

u(x,t) = \tfrac{1}{2}\left[\,F(x - ct) + F(x + ct)\,\right]

$F$ is the initial displacement $f(x)$ extended oddly about both ends with period $2L$ . The fixed-end condition is never imposed separately. The odd extension already flips the sign at each wall. No expansion, no coefficients, no frequencies. The whole solution is the initial shape split in half and carried both ways.

Play with the simulation below.

speed0.35

Turn halves off and you see only the string. Turn it on and the same picture is two rigid copies at half height sliding in opposite directions — no frequencies anywhere. Switch to jump: the corners stay razor sharp forever, because transport never smooths anything. The clock reads t·c/L = 0.00; at 2 the shape is exactly back.

Toggle halves on and off: the single white curve turns out to be two half-height copies sliding past each other. On the jump profile, watch the corners stay razor sharp forever. Advection does not smooth anything.

Bernoulli rewrote the same answer in sines#

Daniel Bernoulli's objection came out of music. One string sounds its fundamental and its overtones at once. So the solution, he argued, has to be a superposition of standing waves.

u(x,t) = \sum_{n=1}^{\infty} b_n \sin\!\frac{n\pi x}{L}\cos\!\frac{n\pi c t}{L}, \qquad b_n = \frac{2}{L}\int_0^L f(x)\,\sin\!\frac{n\pi x}{L}\,dx

$b_n$ is the amplitude of mode $n$ and $n\pi c/L$ is its angular frequency. That the mode frequencies fall in the integer ratio $1:2:3$ had been known since Pythagoras. Paris at the time was arguing about Rameau's theory of harmony, and d'Alembert himself was on the side defending it with Newtonian mechanics. The physical footing of the musical scale came out of this equation.

What Euler actually challenged was the word "function"#

Euler sided with d'Alembert for a different reason. Pluck a string with a finger and the initial shape is a kinked triangle. At the kink, $u_{xx}$ does not exist. d'Alembert did not want to admit such a curve as a solution at all, because for him only a curve expressible as a single analytic formula counted as a function. Euler insisted that any hand-drawn curve is legitimate initial data.

Bernoulli went further still: every curve, he claimed, can be written as an infinite sum of sines. At the time this looked like unfounded optimism. It was not settled until Fourier in 1822. Just as d'Alembert had to live with a potential flow around a cylinder that produces zero drag, his instinct here was half right.

Python, with both representations at the same instant#

Three initial displacements. A box function with jumps, a triangle with only a kink, and a smooth bell. At the same time $tc/L = 0.15$ the d'Alembert solution is compared against an $N$ -term sine series.

import numpy as np
 
L, c, T = 1.0, 1.0, 0.15
X = np.linspace(0.0, L, 2001)
 
def hat_profile(x, a=0.35, b=0.55):
    """Jump: 1 on [a,b], 0 elsewhere - the band a hammer strikes"""
    return np.where((x >= a) & (x <= b), 1.0, 0.0)
 
def kink_profile(x, a=0.45):
    """Kink: continuous, but the slope jumps at x=a - a plucked string"""
    return np.where(x < a, x / a, (L - x) / (L - a))
 
def bell_profile(x, a=0.45, s=0.055):
    """Smooth: infinitely differentiable bell"""
    return np.exp(-((x - a) ** 2) / (2 * s ** 2))
 
def odd_extend(f0, xq):
    """Interpolate on the odd, 2L-periodic extension the fixed ends demand"""
    xs = np.mod(xq, 2 * L)
    sign = np.where(xs > L, -1.0, 1.0)
    xs = np.where(xs > L, 2 * L - xs, xs)
    return sign * np.interp(xs, X, f0)
 
def dalembert_wave(f0, t):
    """d'Alembert solution - two half-height waves carried in opposite directions"""
    return 0.5 * (odd_extend(f0, X - c * t) + odd_extend(f0, X + c * t))
 
def modal_wave(f0, t, n_modes):
    """Bernoulli solution - superposition of n_modes sine modes"""
    n = np.arange(1, n_modes + 1)[:, None]
    k = n * np.pi / L
    b = 2.0 / L * np.trapezoid(f0[None, :] * np.sin(k * X[None, :]), X, axis=1)
    return (b[:, None] * np.sin(k * X[None, :]) * np.cos(k * c * t)).sum(axis=0)
 
def overshoot_pct(u, exact):
    """Maximum excess over the exact peak, as % of the exact amplitude"""
    return 100.0 * (u.max() - exact.max()) / (exact.max() - exact.min())
 
for name, f0 in (("jump  ", hat_profile(X)),
                 ("kink  ", kink_profile(X)),
                 ("smooth", bell_profile(X))):
    exact = dalembert_wave(f0, T)
    print(f"[{name}]  t*c/L = {T}")
    for n_modes in (8, 32, 128, 512):
        u = modal_wave(f0, T, n_modes)
        print(f"   N={n_modes:4d}   max|modal - dAlembert| = {np.abs(u - exact).max():.5f}"
              f"   overshoot = {overshoot_pct(u, exact):+6.2f} %")

[jump  ]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.30480   overshoot = +21.01 %
   N=  32   max|modal - dAlembert| = 0.26454   overshoot = +12.07 %
   N= 128   max|modal - dAlembert| = 0.23847   overshoot =  +9.86 %
   N= 512   max|modal - dAlembert| = 0.18738   overshoot =  +8.94 %
[kink  ]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.02073   overshoot =  -0.23 %
   N=  32   max|modal - dAlembert| = 0.00639   overshoot =  -0.19 %
   N= 128   max|modal - dAlembert| = 0.00157   overshoot =  -0.03 %
   N= 512   max|modal - dAlembert| = 0.00038   overshoot =  -0.01 %
[smooth]  t*c/L = 0.15
   N=   8   max|modal - dAlembert| = 0.04435   overshoot =  -6.25 %
   N=  32   max|modal - dAlembert| = 0.00000   overshoot =  -0.00 %
   N= 128   max|modal - dAlembert| = 0.00000   overshoot =  +0.00 %
   N= 512   max|modal - dAlembert| = 0.00000   overshoot =  -0.00 %

The three blocks tell three different stories. The smooth bell hits double-precision zero by 32 terms. The kinked triangle cuts its error by four every time $N$ goes up by four: first order. Only the jump keeps an overshoot, and it settles near 9% instead of shrinking.

The 9% does not shrink when you add terms#

That number is the Gibbs phenomenon. Wilbraham found it in 1848, Gibbs rediscovered it in 1899 and kept the name. The theoretical value is 8.95% of the jump. The 8.94% printed at 512 terms above is that constant.

The point is that the overshoot narrows without ever getting shorter. Add terms and the ringing region shrinks like $1/N$ . So an integral, or an $L^2$ norm, does converge. Measured in the max norm it does not. That gap between norms is what shows up in practice as a negative density.

modes N24speed0.25

Drag modes N to the right on jump: the red ringing narrows but its height parks near 9% (now 0.00%), and the right-hand curve flattens out. Switch to kink and then smooth with N untouched — the same slider that bought nothing now buys everything, and the error curve turns into a cliff (max error 0.0000).

Push the modes N slider all the way right. On jump the red ringing gets thinner and keeps its height, and the convergence curve on the right flattens along the bottom. Move the same slider on smooth and the curve drops off a cliff. Nothing changed but the initial condition.

The decay rate of the coefficients explains all of it. With a jump, $b_n \sim n^{-1}$ ; with only a kink, $b_n \sim n^{-2}$ ; with a smooth profile, faster than any power. The truncated tail is the error, so the fatter the tail the bigger the scar left at the cut.

1755, when the same man met a nonlinear equation#

Euler wrote fluid motion as a partial differential equation for the first time in 1755: the Euler equations. Unlike the wave equation, here the slope of a characteristic depends on the solution itself. However smooth the initial data, crossing characteristics create a discontinuity in finite time. It is the same structure as why supersonic flow knows nothing about what lies upstream.

So in compressible work, choosing smooth initial data buys you nothing. The shock manufactures its own jump, and from that moment a high-order scheme is back in 1747. von Neumann deliberately smearing shocks in 1950 was aimed at exactly this ringing. TVD limiters and WENO weights drop the order only near a shock, because restoring smoothness locally is the only way to erase the 9%.

Check the smoothness of your initial data first#

When a new scheme rings, look at the initial and boundary data before you suspect the scheme. Was the initial field laid down as cell-wise constants? Was the interface inserted as a jump? Is the inlet profile only $C^0$ in time? If any of those hits, the ringing is not a bug. It is the price of the representation.

The same criterion picks your verification cases. A smooth solution returns the design order cleanly. The moment a jump enters, max-norm convergence disappears and only $L^1$ survives. The three men of 1747 had no vocabulary for that distinction. We call it a choice of norm.