Software 7.5.1 (Coefficients and Parseval). The coefficients of \(f(x) = x\) on \((-\pi , \pi )\) and Parseval’s identity for them; then the coefficients of the plucked string of Example 7.3.1, and their values for \(p = \frac 13\):

>>> import sympy as sp
>>> x, p, h = sp.symbols('x p h', positive=True)
>>> k = sp.symbols('k', integer=True, positive=True)
>>> bk = sp.integrate(x * sp.sin(k*x), (x, -sp.pi, sp.pi)) / sp.pi
>>> bk                                     # f(x) = x on (-pi, pi)
-2*(-1)**k/k
>>> sp.summation(bk**2, (k, 1, sp.oo))     # Parseval: equals (1/pi) int x^2
2*pi**2/3
>>> sp.integrate(x**2, (x, -sp.pi, sp.pi)) / sp.pi
2*pi**2/3
>>> s = sp.sin(k * sp.pi * x)               # plucked at p to height h, L = 1
>>> pk = 2 * (sp.integrate(h*x/p * s, (x, 0, p))
...           + sp.integrate(h*(1 - x)/(1 - p) * s, (x, p, 1)))
>>> sp.factor(sp.simplify(pk))
-2*h*sin(pi*k*p)/(pi**2*k**2*p*(p - 1))
>>> unit = 9 * sp.sqrt(3) / sp.pi**2      # p = 1/3, h = 1: b_k in these units
>>> [sp.simplify(pk.subs({p: sp.Rational(1, 3), h: 1, k: j}) / unit)
...  for j in range(1, 7)]
[1/2, 1/8, 0, -1/32, -1/50, 0]

Reading the output. \(b_k = -2(-1)^k/k\), and \(\sum b_k^2 = 2\pi ^2/3\) equals \(\frac 1\pi \int x^2\), which is \(\sum 1/k^2 = \pi ^2/6\). The plucked-string formula agrees with Example 7.3.1 (the factor \(p - 1\) in the denominator carries the sign). For \(p = \frac 13\) the coefficients in units of \(9\sqrt 3/\pi ^2\) are \(\frac 12, \frac 18, 0, -\frac 1{32}, -\frac 1{50}, 0\): the third and sixth harmonics vanish.

MATLAB/Octave: integral(@(x) x.*sin(k*x), -pi, pi)/pi for each k.

Software 7.5.2 (The Gibbs overshoot). The maximum of the partial sum \(S_N\) of the square wave on \((0, \pi /2]\), and where it occurs, in units of \(1/N\):

>>> import numpy as np
>>> xs = np.linspace(1e-6, np.pi / 2, 200001)          # square wave: sign(x)
>>> def S(N, x):                     # (4/pi) sum, odd k <= N, of sin(kx)/k
...     k = np.arange(1, N + 1, 2)[:, None]
...     return 4 / np.pi * (np.sin(k * x) / k).sum(axis=0)
...
>>> for N in (11, 101, 1001):
...     s = S(N, xs)
...     print(N, round(s.max(), 5), round(float(xs[s.argmax()] * N), 3))
...
11 1.1813 2.88
101 1.17901 3.11
1001 1.17898 3.138
>>> from scipy.special import sici
>>> 2 / np.pi * sici(np.pi)[0]             # the limit of the overshoot
1.178979744472167

Reading the output. The peak value tends to \(1.17898\), the limit \(\frac 2\pi \operatorname {Si}(\pi )\) in the last line; it is \(1.1813\) already with \(11\) terms. The peak sits at \(x \approx \pi /N\) (the products approach \(3.14\)), so the overshoot crowds into the jump without shrinking.

MATLAB/Octave: 2/pi*integral(@(s) sin(s)./s, 0, pi).

Software 7.5.3 (Cooling a rod). A rod of length \(1\) with \(\kappa = 1\), initially at temperature \(1\), with its ends put in ice at \(t = 0\). The sine coefficients of \(g = 1\) are \(4/k\pi \) for odd \(k\). The centre temperature from the series, summed until the coefficients fall below \(10^{-12}\), the last \(k\) used, and the one-term approximation:

>>> def rod(x, t, tol=1e-12):        # u_t = u_xx on (0, 1), u = 0 at ends
...     total, k = 0.0, 1            # u(x, 0) = 1: b_k = 4 / (k pi), k odd
...     while True:
...         c = 4 / (k * np.pi) * np.exp(-(k * np.pi)**2 * t)
...         total += c * np.sin(k * np.pi * x)
...         if c < tol:
...             return total, k
...         k += 2
...
>>> for t in (0.001, 0.01, 0.1, 0.5):
...     u, K = rod(0.5, t)
...     print(t, round(u, 8), K, round(4/np.pi * np.exp(-np.pi**2 * t), 8))
...
0.001 1.0 51 1.26073498
0.01 0.9991861 17 1.15357802
0.1 0.47448746 7 0.47454636
0.5 0.00915699 3 0.00915699

Reading the output. At \(t = 0.001\) the centre has not yet felt the ends (\(1.0\)), and the series needs terms up to \(k = 51\); the single term \(\frac 4\pi e^{-\pi ^2t}\) is useless (\(1.26\)). By \(t = 0.1\) seven terms suffice and one term is within \(10^{-4}\); by \(t = 0.5\) the first mode alone is exact to eight digits. The high modes, which carry the sharp corners of the initial state, are damped first, so for large times the temperature is a pure \(\sin \pi x\) profile decaying at rate \(\pi ^2\kappa /L^2\).

MATLAB/Octave: k = 1:2:51; sum(4./(k*pi).*exp(-(k*pi).^2*t).*sin(k*pi/2)).