Refining the Grid 8x Never Moved 2.4674 — Two Constants Fourier Left Behind in 1807

Refining a grid shrinks the L2 error but never the sawtooth mode's decay, which stays pi^2/4 too slow. Some errors do not respond to refinement at all.

I refined the grid 8x and one number refused to move#

I solved 1D heat conduction on 64 cells. Then 128, 256, 512. The L2 error dropped exactly as expected.

One number printed beside it did not budge. 2.4674. All four grids agreed to four decimal places.

I cut $\Delta x$ by a factor of eight. Why did that error stay put? This post chases that number down. The answer is $\pi^2/4$ , and the trail starts in Paris in 1807.

In 1807 the referees objected to convergence, not physics#

Joseph Fourier submitted a memoir on heat conduction in solids to the Paris Academy in 1807. Its claim: the temperature field inside a solid can be written as a sum of sines and cosines.

Three referees read it — Lagrange, Laplace, Legendre. They rejected it. Not on physical grounds, on mathematical ones. Nobody had shown the series converged.

Lagrange tried to supply the proof himself. He never finished it and died in 1813. The proof arrived much later, from Fourier's own student Dirichlet.

The referees were not being fussy. What they saw was really there. Put a sine series on a function with a jump, and something survives no matter how many terms you add.

On a square pulse of width 0.5, the slider N is how many terms stay in the sum. Raising it packs the spike onto the jump. The height barely moves.

terms N = 41Fejer averagezoom near jump

max overshoot = 0.000%  · peak at x = 0.0000 (jump at x = 0.25)  · distance to jump = -0.2500

Push N up and the spike crowds against the discontinuity. The height sits near 8.95%. Turn on Fejer average and the spike dies, while the step smears.

Stack sines into a square wave and 8.9490% is left over#

Place a square pulse of width 0.5 on a periodic interval. Its Fourier coefficients come out in closed form.

f(x) = \frac{1}{2} + \sum_{m=1}^{\infty} \frac{2}{\pi m}\sin\!\left(\frac{m\pi}{2}\right)\cos\big(2\pi m (x - 0.5)\big)

Here $m$ is the term index and $x$ is position. Only odd $m$ survives. Truncate at $N$ terms and the partial sum shoots above 1 right next to the jump.

This is the Gibbs phenomenon. The limit is 8.9490%, a constant independent of the term count. It is not an error you can buy away with precision.

Why the heat equation is fond of sine series#

Fourier did not pick sines out of taste. In the heat equation

\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}

$\alpha$ is the thermal diffusivity. Substitute $T = e^{ikx}$ and the second derivative returns $-k^2 e^{ikx}$ . So $e^{ikx}$ is an eigenfunction of the Laplacian with eigenvalue $-k^2$ .

Each mode therefore decays on its own as $e^{-\alpha k^2 t}$ , never talking to its neighbors. A partial differential equation splits into one ordinary differential equation per mode.

That $k^2$ inside the exponent sets the decay rate. Halve the wavelength and the decay is four times faster. Short waves get erased first. Everything diffusion does is packed into that one spot.

What the grid hands back instead of $-k^2$ #

Now run the same calculation on a grid. Replace the second derivative with a second-order central difference and substitute $T_j = e^{ikj\Delta x}$ .

\frac{T_{j+1} - 2T_j + T_{j-1}}{\Delta x^2} = -\frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right) e^{ikj\Delta x}

Here $j$ is the node index. The discrete operator also has $e^{ikj\Delta x}$ as an eigenfunction, which means Fourier decomposition survives the discretization intact.

What changed is the eigenvalue. Not $-k^2$ but $-k'^2$ , with

k'^2 = \frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right)

This $k'$ is the modified wavenumber. It is the scheme confessing which $k$ it thinks it is being given.

For small $\theta = k\Delta x$ the effective curve sits on the $\theta^2$ parabola. Toward the right edge the two split. The same wave dies later on the grid than in the exact exponent.

Where $\pi^2/4$ shows up at the Nyquist limit#

The shortest wavelength a grid can carry is $2\Delta x$ : up one cell, down the next, a sawtooth. There $\theta = k\Delta x = \pi$ .

Feed that into the formula above:

\frac{k'^2}{k^2} = \frac{4\sin^2(\pi/2)}{\pi^2} = \frac{4}{\pi^2} = 0.4053

The reciprocal is $\pi^2/4 = 2.4674$ — the number that would not move.

The important detail is that $\Delta x$ cancelled. The ratio does not depend on the mesh spacing. However tightly you refine, that grid's sawtooth always dies 2.4674 times too slowly.

Raising the order barely helps. Here it is on real stencils.

import numpy as np
 
STENCILS = {
    2: [1.0, -2.0, 1.0],
    4: [-1/12, 4/3, -5/2, 4/3, -1/12],
    6: [1/90, -3/20, 3/2, -49/18, 3/2, -3/20, 1/90],
}
 
def modified_k2(order, theta):
    """Effective squared wavenumber of the central 2nd derivative (dx=1). Exact value is theta**2."""
    c = STENCILS[order]
    h = len(c) // 2
    sym = sum(c[h + j] * np.cos(j * theta) for j in range(-h, h + 1))
    return -sym
 
theta = np.array([np.pi/8, np.pi/4, np.pi/2, 3*np.pi/4, np.pi])
print("theta/pi   exact k^2    order2    order4    order6    (ratio)")
for th in theta:
    ex = th**2
    row = [modified_k2(o, th) for o in (2, 4, 6)]
    print(f"{th/np.pi:6.3f}  {ex:10.5f}  " + "  ".join(f"{v:7.4f}" for v in row)
          + "   " + " ".join(f"{v/ex:6.4f}" for v in row))
 
print()
for o in (2, 4, 6):
    r = modified_k2(o, np.pi) / np.pi**2
    print(f"order {o}: Nyquist ratio = {r:.4f}   decay is {1/r:.3f}x slower")

theta/pi   exact k^2    order2    order4    order6    (ratio)
 0.125     0.15421   0.1522   0.1542   0.1542   0.9872 0.9997 1.0000
 0.250     0.61685   0.5858   0.6144   0.6166   0.9496 0.9960 0.9996
 0.500     2.46740   2.0000   2.3333   2.4222   0.8106 0.9457 0.9817
 0.750     5.55165   3.4142   4.3856   4.8278   0.6150 0.7900 0.8696
 1.000     9.86960   4.0000   5.3333   6.0444   0.4053 0.5404 0.6124
 
order 2: Nyquist ratio = 0.4053   decay is 2.467x slower
order 4: Nyquist ratio = 0.5404   decay is 1.851x slower
order 6: Nyquist ratio = 0.6124   decay is 1.633x slower

The dial repeats that table. Order lifts the low-wavenumber side onto $\theta^2$ . The sawtooth ratio barely moves.

theta/pi = 1.00

k'^2 / k^2 = 0.4053  · grid decays 2.467x slower  · a*t/dx^2 = 0.000

Push $\theta$ out to $\pi$ and the grid copy decays later than the exact wave. Sixth order still sits near a Nyquist ratio of 0.6124, which is 1.633 times too slow.

Waves spanning four cells or more ( $\theta/\pi \le 0.25$ ) become nearly perfect once the order goes up. The sawtooth does not. Even sixth order reads 0.6124.

The same square wave on four grids#

Now actually solve it. To keep time-integration error out of the picture, the time direction is integrated analytically, leaving only the spatial discretization.

import numpy as np
 
L, ALPHA, TEND = 1.0, 1.0, 2.0e-3
 
def tophat_coeffs(m):
    """Complex Fourier coefficients of the top hat 0.25 <= x < 0.75 on the period [0,1)."""
    m = np.asarray(m, dtype=float)
    out = np.zeros(m.shape, dtype=complex)
    out[m == 0] = 0.5
    mm = m[m != 0]
    out[m != 0] = (np.exp(-2j*np.pi*mm*0.75) - np.exp(-2j*np.pi*mm*0.25)) / (-2j*np.pi*mm)
    return out
 
def gibbs_overshoot(nterm):
    """Maximum overshoot (%) of the 2N+1 term partial sum. Sampled densely only near the spike."""
    x = np.linspace(0.25, 0.25 + 4.0/nterm, 40001)
    m = np.arange(-nterm, nterm + 1)
    c = tophat_coeffs(m)
    s = (c[None, :] * np.exp(2j*np.pi*m[None, :]*x[:, None])).sum(axis=1).real
    return (s.max() - 1.0) * 100.0
 
def heat_error(ngrid, smooth):
    """Heat equation integrated exactly in time; only the spatial discretization differs."""
    dx = L / ngrid
    x = np.arange(ngrid) * dx
    m = np.fft.fftfreq(ngrid, d=1.0/ngrid)
    k = 2*np.pi*m
    kmod2 = (4.0/dx**2) * np.sin(k*dx/2.0)**2          # k^2 the grid actually delivers
    mref = np.arange(-8192, 8193)
    if smooth:
        t0 = np.cos(2*np.pi*x) + 0.3*np.cos(6*np.pi*x)
        ref = (np.cos(2*np.pi*x)*np.exp(-ALPHA*(2*np.pi)**2*TEND)
               + 0.3*np.cos(6*np.pi*x)*np.exp(-ALPHA*(6*np.pi)**2*TEND))
    else:
        t0 = ((x >= 0.25) & (x < 0.75)).astype(float)
        c = tophat_coeffs(mref)*np.exp(-ALPHA*(2*np.pi*mref)**2*TEND)
        ref = (c[None, :]*np.exp(2j*np.pi*mref[None, :]*x[:, None])).sum(axis=1).real
    num = np.fft.ifft(np.fft.fft(t0)*np.exp(-ALPHA*kmod2*TEND)).real
    nyq = ngrid//2
    return np.sqrt(np.mean((num-ref)**2)), 1.0/(ALPHA*k[nyq]**2), 1.0/(ALPHA*kmod2[nyq])
 
print("N_terms   Gibbs overshoot[%]")
for n in (16, 32, 64, 128, 256, 512):
    print(f"{n:6d}    {gibbs_overshoot(n):9.5f}")
 
for smooth in (True, False):
    print(f"\ninitial condition: {'smooth cosine' if smooth else 'square wave (discontinuous)'}")
    print("grid      L2 error     order     sawtooth decay time exact/grid    ratio")
    prev = None
    for ng in (64, 128, 256, 512):
        l2, t_ex, t_gr = heat_error(ng, smooth)
        order = "  -  " if prev is None else f"{np.log2(prev/l2):5.2f}"
        print(f"{ng:4d}   {l2:.4e}    {order}     {t_ex:.3e} / {t_gr:.3e}   {t_gr/t_ex:.4f}")
        prev = l2

N_terms   Gibbs overshoot[%]
    16      9.01421
    32      8.96527
    64      8.95306
   128      8.95000
   256      8.94924
   512      8.94905
 
initial condition: smooth cosine
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   5.3683e-04      -       2.474e-05 / 6.104e-05   2.4674
 128   1.3424e-04     2.00     6.184e-06 / 1.526e-05   2.4674
 256   3.3562e-05     2.00     1.546e-06 / 3.815e-06   2.4674
 512   8.3907e-06     2.00     3.865e-07 / 9.537e-07   2.4674
 
initial condition: square wave (discontinuous)
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   2.3383e-02      -       2.474e-05 / 6.104e-05   2.4674
 128   1.1673e-02     1.00     6.184e-06 / 1.526e-05   2.4674
 256   5.8342e-03     1.00     1.546e-06 / 3.815e-06   2.4674
 512   2.9168e-03     1.00     3.865e-07 / 9.537e-07   2.4674

Errors that shrink with refinement, and errors that do not#

With the smooth initial condition the convergence order is exactly 2.00. The scheme is second order, as advertised.

Hand the same scheme a square wave and it drops to 1.00. A discontinuity makes the Fourier coefficients decay only as $1/m$ , and that slow tail returns as aliasing. Convergence order is set by the smoothness of the solution, not by the order of the scheme. Killing the ringing next to a step is the same job as building a second-order TVD scheme with a flux limiter.

The ratio in the last column is 2.4674 on every one of the eight rows. The L2 error fell by a factor of eight; this value is identical past the decimal point.

Do not blur those two facts. L2 falls because the sawtooth mode's amplitude gets smaller, not because the scheme got better at handling sawtooths. Refining a grid does not solve the problem better; it moves the problem.

On 64 cells the smooth cosine has L2 equal to 5.3683e-04. Eight times finer, 512 cells, that is 8.3907e-06. The ratio is about 64, and the log still says order 2.00. The square wave on the same grids falls from 2.3383e-02 to 2.9168e-03, a factor of eight, order 1.00. Sawtooth decay time at 64 cells is 2.474e-05 exact against 6.104e-05 on the grid. At 512 cells the pair is 3.865e-07 and 9.537e-07. The ratio stays 2.4674 through the fourth decimal. Exact $k^2$ at Nyquist is 9.86960. The second-order stencil returns 4.0000, and $4/\pi^2 = 0.4053$ . Once $\Delta x$ cancels, refinement has nothing left to change.

Gibbs behaves the same way. At 16 terms the overshoot is 9.01421%. At 512 terms it is 8.94905%. The gap is about 0.065 percentage points. On the $\theta/\pi = 0.125$ row, order 2 already reads 0.9872, so sixteen cells per wavelength sits inside 1%. Sixth order at Nyquist returns 6.0444, ratio 0.6124, decay 1.633 times too slow. Buying order buys digits at low wavenumber. It does not buy the sawtooth.

Divide the square-wave L2 at 64 cells, 2.3383e-02, by the smooth L2, 5.3683e-04. About 44. At 512 cells the pair is 2.9168e-03 against 8.3907e-06, about 348. Refinement makes the discontinuous case relatively worse, because orders 1 and 2 separate by one power of the grid ratio. 512/64 = 8, so the error ratio should widen by about 8. 44 times 8 is 352, and the table says 348.

Why sixth order still reads 0.6124, and what filters are for#

Raising the order only pays at low wavenumbers. On the $\theta/\pi = 0.125$ row sixth order gives 1.0000; at Nyquist it gives 0.6124. No finite stencil handles the sawtooth correctly, because differentiation of a $2\Delta x$ wave is genuinely ambiguous on the grid.

So practice splits two ways. One is to keep the sawtooth from forming at all. Staggered grids and Rhie–Chow interpolation do that work. Tightening only the bad cells still leaves that cell's sawtooth ratio at $\pi^2/4$ . A quadtree that grows the mesh on its own does not erase the ratio either.

The other is to erase the sawtooth explicitly once it exists. Spectral filters, artificial viscosity and WENO weights all belong here. The Fejer average in the viz above is the simplest member of the family, and it shows the price in the same picture: the Gibbs spike goes, and the step smears.

Neither route is free. What you choose is not whether the error exists, but where you put it.

The blank Dirichlet filled in, in 1829#

Lagrange's objection was not a formality. He doubted what a series does when it meets a jump, and 8.9490% really was sitting there.

Dirichlet proved conditional convergence in 1829 and filled the blank. The series converges at each point — just not uniformly. The Gibbs spike was never a counterexample; it was what "not uniformly" looks like.

Every grid-convergence study repeats that question. Is this error the kind that shrinks with refinement, or the kind where $\Delta x$ cancels and $\pi^2/4$ stays? A single row of convergence order answers only the first.

If a clean 2.00 has ever put you at ease, plotting the modified wavenumber once is worth the trouble. What a scheme confesses about the sawtooth is written nowhere else.