Convergence of Fourier series and Gibbs phenomenon

Lesson 5 computed Fourier series without asking whether they actually add up to the function. That question is subtle. There are continuous functions whose Fourier series diverge at some points, a discovery of Paul du Bois-Reymond in 1873, and the full story occupied analysts for more than a century. Fortunately, for the functions that arise in PDE problems, piecewise smooth functions with finitely many jumps and corners, the answer is clean. This lesson states and proves Dirichlet's theorem on pointwise convergence, explains when convergence is uniform and when it holds only in the mean-square sense, relates smoothness to the decay of coefficients, and analyses the Gibbs phenomenon, the stubborn overshoot of partial sums next to a jump.

Piecewise smooth functions

A function \(f\) on \([a, b]\) is piecewise continuous if it is continuous except at finitely many points, where it has finite one-sided limits \(f(x^+)\) and \(f(x^-)\). It is piecewise smooth if both \(f\) and \(f'\) are piecewise continuous. The square wave, the triangle wave and the sawtooth \(f(x) = x\) are all piecewise smooth; \(\sqrt{|x|}\) is not, because its derivative is unbounded at 0, and \(x\sin(1/x)\) is not, because its derivative oscillates without limit. Throughout, \(f\) on \([-L, L]\) is extended to be \(2L\)-periodic, which may create jumps at \(\pm L\) even when \(f\) is continuous on the interval.

Dirichlet's theorem

In 1829 Peter Gustav Lejeune Dirichlet published the first rigorous convergence theorem, settling the question for the functions Fourier had in mind.

for every \(x\), when \(f\) is piecewise smooth and \(2L\)-periodic. At points of continuity the limit is \(f(x)\); at a jump it is the midpoint of the jump. So the square wave series converges to \(\pm1\) away from the jumps and to 0 at them, exactly as the figure in lesson 5 showed. Values of \(f\) at isolated points have no effect on the coefficients, so the series cannot tell what value \(f\) was given at a jump; the average is the natural compromise.

Proof by the Dirichlet kernel

Take \(L = \pi\) for simplicity. Substituting the coefficient formulas into \(S_N\) and using \(\cos n(x - t) = \cos nx\cos nt + \sin nx\sin nt\) gives

after the substitution \(s = t - x\) and periodicity. The Dirichlet kernel has a closed form, found by multiplying by \(2\sin(s/2)\) and telescoping with \(2\sin(s/2)\cos ns = \sin(n + \frac12)s - \sin(n - \frac12)s\):

The kernel integrates to 1 over \([-\pi, \pi]\) (integrate the cosine sum term by term), and it is even, so it integrates to \(\frac12\) over each half \([0, \pi]\) and \([-\pi, 0]\). As \(N\) grows, it develops a tall central spike of height \((2N + 1)/(2\pi)\) and oscillates ever faster elsewhere. Splitting the integral into halves and subtracting the target,

In the first integral write \(D_N(s) = \sin((N + \frac12)s)\cdot g(s)\) with \(g(s) = \frac{f(x + s) - f(x^+)}{2\pi\sin(s/2)}\). Because \(f\) is piecewise smooth, the quotient \((f(x + s) - f(x^+))/s\) stays bounded as \(s \to 0^+\) (it tends to the one-sided derivative), and \(s/\sin(s/2)\) is bounded too. So \(g\) is piecewise continuous on \([0, \pi]\). The integral is then the Fourier sine coefficient of an integrable function at frequency \(N + \frac12\), which tends to zero by the Riemann–Lebesgue lemma: for any integrable \(g\), \(\int g(s)\sin(\lambda s)\,ds \to 0\) as \(\lambda \to \infty\), because rapid oscillation cancels the contributions of nearby points. The same argument handles the second integral, completing the proof.

The proof shows where smoothness is used. All that matters is that \(f\) is integrable and has one-sided derivatives at the point in question, a local condition. This is Riemann's localization principle: convergence of a Fourier series at a point depends only on the behavior of \(f\) near that point, even though every coefficient depends on \(f\) everywhere.

Uniform convergence, smoothness and decay

Pointwise convergence can be uneven, and for PDE work we often need more. A series converges uniformly if \(\max_x |f(x) - S_N(x)| \to 0\). The partial sums are continuous, and a uniform limit of continuous functions is continuous, so a function with a jump can never be the uniform limit of its Fourier series. The converse holds under mild conditions: if \(f\) is continuous, \(2L\)-periodic (so no jump at \(\pm L\)) and piecewise smooth, its Fourier series converges absolutely and uniformly to \(f\). The proof uses the coefficients of \(f'\). Integration by parts gives \(a_n(f') = (n\pi/L)\,b_n(f)\) and \(b_n(f') = -(n\pi/L)\,a_n(f)\), with no boundary terms because \(f\) is continuous and periodic. By the Cauchy–Schwarz inequality and Bessel's inequality for \(f'\),

and the Weierstrass M-test gives uniform convergence; Dirichlet's theorem identifies the limit as \(f\). The same integration by parts shows the general pattern. Each derivative that is continuous and periodic adds one power of \(n\) to the decay of the coefficients. A function with jumps has coefficients of order \(1/n\); a continuous function with corners, order \(1/n^2\); a function with continuous first derivative and jumps in the second, \(1/n^3\); and an infinitely smooth periodic function has coefficients that decay faster than any power of \(n\), often exponentially. Conversely, slow decay of the coefficients signals roughness in the function.

Term-by-term operations follow the same rules. A Fourier series can always be integrated term by term, and the result converges uniformly, because integration improves decay. It can be differentiated term by term only when \(f\) is continuous, periodic and piecewise smooth with \(f'\) piecewise smooth, because differentiation worsens decay. Differentiating the square wave series term by term gives \(\frac4\pi(\cos x + \cos 3x + \cdots)\), which diverges everywhere.

Mean-square convergence

For merely square-integrable functions there is always convergence in the mean-square sense: \(\int_{-L}^L (f - S_N)^2\,dx \to 0\). This is equivalent to Parseval's identity, and it says that the trigonometric system is complete, meaning that no nonzero function is orthogonal to all of its members. Mean-square convergence does not imply convergence at any particular point, but it is exactly the kind of convergence that energy arguments in PDE require. The completeness of the trigonometric system, together with the theory of the Lebesgue integral, makes the space of square-integrable functions a Hilbert space with the trigonometric functions as an orthonormal basis, a structure that underlies quantum mechanics and much of modern analysis. Lennart Carleson proved in 1966 that the Fourier series of every square-integrable function converges at almost every point, a result that had been an open problem for over a century.

The Gibbs phenomenon

Near a jump, the partial sums of a Fourier series overshoot. For the square wave, the maximum of \(S_N\) just to the right of 0 does not approach 1 as \(N\) grows; it approaches about 1.179. The overshoot moves closer to the jump and gets narrower, but its height stays the same. The effect was first described by Henry Wilbraham in 1848, then forgotten. It became famous in 1898, when Albert Michelson, who had built a mechanical harmonic analyzer that summed 80 terms, wrote to Nature about apparent anomalies in the sum of the sawtooth series, and Josiah Willard Gibbs explained the limiting curve in letters of 1898 and 1899. Maxime Bôcher gave the full analysis in 1906 and named it after Gibbs.

The height can be computed exactly. For the square wave, \(S_N\) with harmonics up to \(N = 2M - 1\) satisfies \(S_N'(x) = \frac4\pi\sum_{k=1}^{M}\cos(2k - 1)x = \frac{2\sin(2Mx)}{\pi\sin x}\), whose first positive zero is \(x = \pi/(2M)\). At that point the sum is a Riemann sum:

Here \(\mathrm{Si}\) is the sine integral. The jump in the square wave has size 2, so the overshoot above the value 1 is 0.17898, about 8.95 hundredths of the jump. The same fraction appears at every jump of every piecewise smooth function, because near a jump the function looks like a scaled square wave plus a continuous function whose series converges uniformly.

The Gibbs phenomenon is not a failure of Dirichlet's theorem. At each fixed \(x > 0\), the peak eventually moves past \(x\) and \(S_N(x) \to 1\). What fails is uniform convergence: the maximum error near the jump does not shrink. The phenomenon matters in practice. Truncated Fourier representations of images produce ringing artifacts near sharp edges; spectral methods for PDEs with shocks produce spurious oscillations; and filters in signal processing that cut off sharply in frequency produce overshoot in time.

Taming the overshoot

The remedy is to average or to taper the coefficients. Lipót Fejér showed in 1900 that the arithmetic means of the partial sums, \(\sigma_N = (S_0 + S_1 + \cdots + S_{N-1})/N\), converge uniformly to every continuous periodic function and never overshoot: their kernel is non-negative. In terms of coefficients, the Fejér mean multiplies the \(n\)-th harmonic by \(1 - n/N\). Cornelius Lanczos proposed the sigma factors \(\mathrm{sinc}(n\pi/N) = \sin(n\pi/N)/(n\pi/N)\), which reduce the overshoot to about 1.2 hundredths of the jump while keeping the edge sharper. The price of any such filter is a slightly smeared jump.

Measuring Gibbs numerically

The program below locates the first peak of the square wave partial sums for increasing numbers of terms, compares its height with \(\frac2\pi\mathrm{Si}(\pi)\) computed by Simpson's rule, measures the maximum error near and away from the jump, and shows the effect of Fejér and Lanczos smoothing.

import math

def S(x, N, w=lambda n: 1.0):
    # square wave partial sum with optional coefficient weights w(n)
    return 4 / math.pi * sum(w(n) * math.sin(n * x) / n for n in range(1, N + 1, 2))

def simpson(g, a, b, n=1000):
    h = (b - a) / n
    return h / 3 * (g(a) + g(b) + sum((4 if i % 2 else 2) * g(a + i * h) for i in range(1, n)))

Si_pi = simpson(lambda t: math.sin(t) / t if t else 1.0, 0, math.pi)
print(f"(2/pi) Si(pi) = {2 / math.pi * Si_pi:.6f}")

for N in (11, 41, 161, 641):
    xp = math.pi / (N + 1)                    # first peak of S_N
    peak = S(xp, N)
    away = max(abs(1 - S(x, N)) for x in [0.5 + 0.01 * k for k in range(200)])
    print(f"N = {N:3d}: peak at x = {xp:.5f}, height {peak:.5f}, max error on [0.5, 2.5] = {away:.5f}")

N = 161
fejer = lambda n: 1 - n / (N + 1)
lanczos = lambda n: math.sin(n * math.pi / (N + 1)) / (n * math.pi / (N + 1))
grid = [k * 0.0005 for k in range(1, 2001)]
for name, w in (("plain  ", lambda n: 1.0), ("Fejer  ", fejer), ("Lanczos", lanczos)):
    top = max(S(x, N, w) for x in grid)
    rise = next(x for x in grid if S(x, N, w) > 0.9)
    print(f"{name}: max over (0, 1] = {top:.5f}, reaches 0.9 at x = {rise:.4f}")

# convergence at the jump itself and decay of coefficients versus smoothness
print("S_N(0) =", S(0.0, 641))
for name, c in (("square, 1/n", lambda n: 4 / (math.pi * n)), ("triangle, 1/n^2", lambda n: 4 / (math.pi * n * n))):
    print(name, [f"{c(n):.2e}" for n in (1, 11, 101, 1001)])

(2/pi) Si(pi) = 1.178980
N =  11: peak at x = 0.26180, height 1.18130, max error on [0.5, 2.5] = 0.10184
N =  41: peak at x = 0.07480, height 1.17917, max error on [0.5, 2.5] = 0.02985
N = 161: peak at x = 0.01939, height 1.17899, max error on [0.5, 2.5] = 0.00686
N = 641: peak at x = 0.00489, height 1.17898, max error on [0.5, 2.5] = 0.00176
plain  : max over (0, 1] = 1.17896, reaches 0.9 at x = 0.0105
Fejer  : max over (0, 1] = 0.99535, reaches 0.9 at x = 0.0330
Lanczos: max over (0, 1] = 1.02371, reaches 0.9 at x = 0.0195
S_N(0) = 0.0
square, 1/n ['1.27e+00', '1.16e-01', '1.26e-02', '1.27e-03']
triangle, 1/n^2 ['1.27e+00', '1.05e-02', '1.25e-04', '1.27e-06']

The peak height settles at 1.1790 to four decimals while the maximum error away from the jump shrinks steadily, which is pointwise but not uniform convergence. The Fejér mean never exceeds 1 but rises more slowly; the Lanczos factors leave a small overshoot and a sharper edge. At the jump itself the partial sum is exactly 0, the average of the one-sided limits.