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

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.

Toggle halves off and the white curve splits into two half-height travelers. Switch to jump and the corners stay sharp.

The string is the sum of those two halves. The profile that keeps its corners is d'Alembert's solution. The modal sum leaves the 9% for later.

The 9% ringing 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$ 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.

$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. 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.

$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. A linear wave only carries its shape. Once characteristics cross, a nonlinear problem grows its own shock even from smooth data. d'Alembert's instinct was only half right here.

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 %

Sweep modes N from 8 to 512. On jump the overshoot stalls near 9% while the red ringing thins out. On smooth the error is gone by 32. On kink the error drops by about four each time N is multiplied by four.

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 jump's max error falls from 0.30480 at $N = 8$ to 0.18738 at $N = 512$. Sixty-four times more terms buy only a factor of about 1.6. The overshoot drops from +21.01% to +8.94%, within 0.01 points of the theoretical 8.95%. Over the same range the kink prints 0.02073, 0.00639, 0.00157, 0.00038. The ratios are about 3.24, 4.07 and 4.13, close to the factor-of-four drop expected when $N$ grows by four. The bell has error 0.04435 and overshoot −6.25% at $N = 8$, then 0.00000 at $N = 32$. At $tc/L = 0.15$ the traveling waves have moved only 15% of the length. The scar from the truncated series is already sitting on the corner.

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 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. 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%.

Max-norm order disappears the moment the data jumps#

When a new scheme rings, do not open the limiter first. Check whether the initial field is cell-wise constant, whether the interface was inserted as a jump, and whether the inlet is only $C^0$ in time. Data whose overshoot is still +8.94% at 512 terms is paying for the representation.

Verification cases use the same cut. The bell returns the design order outright: error 0.00000 by $N = 32$. The moment a jump enters, max-norm convergence disappears and only $L^1$ survives. The three men of 1747 had no name for that split. The name is the choice of norm.