The Gibbs overshoot never shrinks

Take the square wave that is $+1$ on $(0,\pi )$ and $-1$ on $(-\pi ,0)$. Its Fourier series is $\frac{4}{\pi }(\sin x+\frac{1}{3}\sin 3x+\frac{1}{5}\sin 5x+\dots )$: odd harmonics only, falling off like $1/k$. Sum the first five and you get a wobbly approximation. Sum eighty and the wobble has all but gone — except right next to the jump, where a spike of the same height as ever has moved in toward $x=0$ and thinned.

Where the height comes from

The sum of the first $m$ odd harmonics has its first peak at exactly $x=\pi /2m$, where the derivative $\frac{2}{\pi }\sin (2mx)/\sin x$ first vanishes. Evaluate the sum there and it is a Riemann sum, with $m$ panels, for $\frac{2}{\pi }{\int }_0^{\pi }\frac{\sin t}{t}dt=1.17898$. So the height of the first peak tends to a number that has nothing to do with $m$: $(2/\pi )\mathrm{Si}(\pi )$, against a wave of height 1. The overshoot is $0.179$, which is $8.95\%$ of the jump from $-1$ to $1$.

Why more terms cannot help

The peak sits at $\pi /2m$ and is about that wide, so as $m$ grows it carries less and less area. In the mean-square sense the series converges: the integral of the squared error goes to zero. But at the peak itself the error stays $0.179$ for every $m$. That is the difference between convergence in $L^2$ and uniform convergence, and a square wave, being discontinuous, has the first without the second. A smoothing filter (Fejér, Lanczos) trades the spike for a slower approach to the wave; nothing that keeps the partial sums as they are removes it.

What the program checks

The program checks three things. The peak of a thousand harmonics, summed directly, matches the integral computed by Simpson's rule with 2000 panels to a millionth. The peaks of ten and of a hundred harmonics are both above $1.17$, so adding terms does not lower the peak. And the overshoot is $0.0895$ of the jump, to four places.