Software 4.6.1 (D’Alembert and Duhamel). A symbolic check of (4.2) for arbitrary \(f\) and \(g\), and the Duhamel integral of Theorem 4.4.1 for the source \(F = xt\):

>>> import sympy as sp
>>> x, t, s, y, c = sp.symbols('x t s y c', positive=True)
>>> f, g = sp.Function('f'), sp.Function('g')
>>> u = ((f(x + c*t) + f(x - c*t)) / 2
...      + sp.Integral(g(y), (y, x - c*t, x + c*t)) / (2*c))
>>> sp.simplify(u.diff(t, 2) - c**2 * u.diff(x, 2))          # d'Alembert
0
>>> u.subs(t, 0).doit(), sp.simplify(u.diff(t).subs(t, 0).doit())
(f(x), g(x))
>>> F = lambda y, s: y * s                 # forcing in u_tt = c^2 u_xx + F
>>> inner = sp.integrate(F(y, s), (y, x - c*(t - s), x + c*(t - s)))
>>> w = sp.integrate(inner, (s, 0, t)) / (2*c)   # Duhamel: over the triangle
>>> sp.factor(w), sp.simplify(w.diff(t, 2) - c**2 * w.diff(x, 2))
(t**3*x/6, t*x)

Reading the output. D’Alembert’s formula satisfies the equation identically and takes the right initial values. The Duhamel integral over the triangle gives \(u = xt^3/6\), and \(u_{tt} - c^2u_{xx} = xt\) returns the source.

MATLAB/Octave: with the Symbolic Math Toolbox, int(int(y*s, y, x-c*(t-s), x+c*(t-s)), s, 0, t)/(2*c).

Software 4.6.2 (Energy of a plucked string). A string of length \(1\), with \(c = 1\) and fixed ends, plucked at \(x = 0.3\) to a height of \(0.01\) and released, is solved by the leapfrog scheme \(u_j^{n+1} = 2u_j^n - u_j^{n-1} + r^2(u_{j+1}^n - 2u_j^n + u_{j-1}^n)\) with \(r = c\,\Delta t/\Delta x = \frac 12\). The energy is computed in the form that the scheme conserves exactly, \(\frac 12\sum \bigl ((u_j^{n+1} - u_j^n)/\Delta t\bigr )^2\Delta x + \frac 12\sum (u_{j+1}^{n+1} - u_j^{n+1})(u_{j+1}^n - u_j^n)/\Delta x\):

>>> import numpy as np
>>> def leapfrog(u0, r, steps, ends='fixed'):   # u_tt = u_xx, r = c dt / dx
...     uo = u0.copy(); un = u0.copy()
...     un[1:-1] += 0.5 * r**2 * (u0[2:] - 2*u0[1:-1] + u0[:-2])
...     for _ in range(steps - 1):
...         nxt = np.empty_like(un)
...         nxt[1:-1] = (2*un[1:-1] - uo[1:-1]
...                      + r**2 * (un[2:] - 2*un[1:-1] + un[:-2]))
...         nxt[0] = 0 if ends == 'fixed' else nxt[1]      # u = 0 or u_x = 0
...         nxt[-1] = 0
...         uo, un = un, nxt
...     return uo, un
...
>>> x = np.linspace(0, 1, 401); dx = x[1]                  # string, L = c = 1
>>> u0 = 0.01 * np.where(x < 0.3, x / 0.3, (1 - x) / 0.7)  # plucked at 0.3
>>> 0.5 * np.sum(np.diff(u0)**2) / dx      # energy = (1/2) int u_x^2 at rest
0.00023809523809523812
>>> for steps in (400, 800, 1600):         # t = 0.5, 1, 2 with r = 1/2
...     uo, un = leapfrog(u0, 0.5, steps); dt = 0.5 * dx
...     E = 0.5 * np.sum(((un - uo) / dt)**2) * dx \
...         + 0.5 * np.sum(np.diff(un) * np.diff(uo)) / dx
...     print(steps * dt, f"{E:.6e}", round(float(np.abs(un - u0).max()), 6))
...
0.5 2.379181e-04 0.011906
1.0 2.379181e-04 0.014325
2.0 2.379181e-04 0.000128

Reading the output. The exact energy is \(\frac 12\int f'^2 = \frac 12(0.01)^2(\frac 1{0.3} + \frac 1{0.7}) = 2.381 \times 10^{-4}\). The computed energy, \(2.379181 \times 10^{-4}\), is the same to all seven digits at \(t = 0.5\), \(1\) and \(2\); it differs from the exact value by \(0.07\%\), an error made by the starting step. The last column is the largest distance from the initial shape: large at \(t = 0.5\) and \(t = 1\), and \(1.3 \times 10^{-4}\), about \(1\%\) of the amplitude, at \(t = 2 = 2L/c\), when the string has returned to where it started.

MATLAB/Octave: the same loop with u(2:end-1) indexing.

Software 4.6.3 (Reflection at fixed and free ends). A pulse \(e^{-50(x - 2)^2}\), released from rest on the half-line, with a fixed end (\(u = 0\)) and with a free end (\(u_x = 0\)) at \(x = 0\); the solution at \(t = 3\):

>>> x = np.linspace(0, 6, 1201); dx = x[1]  # half-line, cut at 6
>>> pulse = lambda z: np.exp(-50 * (z - 2)**2)  # starts at rest at 2
>>> for ends in ('fixed', 'free'):
...     _, u = leapfrog(pulse(x), 0.5, 1200, ends)  # t = 1200 dx / 2 = 3
...     print(ends, np.round(np.interp([0.5, 1.0, 1.5, 5.0], x, u), 4))
...
fixed [-0.  -0.5 -0.   0.5]
free [0.     0.4995 0.     0.5   ]

Reading the output. The pulse splits into two halves of height \(\frac 12\). The right-moving half is at \(x = 2 + 3 = 5\). The left-moving half reached the end at \(t = 2\) and has come back to \(x = 1\), upside down (\(-0.5\)) at the fixed end and the right way up (\(0.4995\)) at the free end, as the odd and even extensions predict.

MATLAB/Octave: the same scheme, with u(1) = u(2) for the free end.