The Stretched String Problem

We’ll start our discussion of numerical methods by considering the problem of finding normal-mode eigenfrequencies and eigenfunctions for waves on a stretched string clamped at both ends. Let the string have mass per unit length \(\rho\) and tension \(T\); then, the wave equation describing the transverse string displacement \(y(x,t)\) at spatial position \(x\) and time \(t\) is

\[\npderiv{y}{x}{2} = \frac{1}{c^{2}} \npderiv{y}{t}{2},\]

with \(c \equiv (T/\rho)^{1/2}\). If the string is clamped at \(x=0\) and \(x=L\), then the wave equation together with the boundary conditions

\[y(0,t) = 0 \qquad y(L,t) = 0\]

comprises a two-point boundary value problem (BVP).

Analytic Solution

The stretched-string BVP is straightforward to solve analytically. General solutions of the wave equation take the form of traveling waves,

\[y(x,t) = A \exp [\ii (k x - \sigma t) ],\]

where \(A\) an arbitrary constant, and the frequency \(\sigma\) and wavenumber \(k\) are linked by the dispersion relation

\[\sigma^{2} = c^{2} k^{2}.\]

The phase velocity of these waves is \(\sigma/k = \pm c\).

To satisfy the boundary condition at \(x=0\), we combine traveling-wave solutions with opposite-sign wavenumbers

\[y(x,t) = A \exp [\ii (k x - \sigma t) ] - A \exp [\ii (- k x - \sigma t) ] = B \sin(k x) \exp ( - \ii \sigma t),\]

where \(B = 2A\). For the boundary condition at \(x=L\) to be satisfied simultaneously,

\[\sin(k L) = 0,\]

and so

\[k L = n \pi\]

where \(n\) is a non-zero integer (we exclude \(n=0\) because it corresponds to the trivial solution \(y(x,t)=0\)). Combining this with the dispersion relation, we find that the normal-mode eigenfrequencies of the stretched-string BVP are

(1)\[\sigma = n \frac{\pi c}{L},\]

and the corresponding eigenfunctions are

(2)\[y_{n}(x,t) = B \sin \left( \frac{n \pi x}{L} \right) \exp ( - \ii \sigma t).\]

The index \(n\) uniquely labels the modes, and \(y_{n}(x,t)\) exhibits \(n-1\) nodes in the open interval \(x \in (0,L)\).

Separation

Now let’s see how we might go about solving the stretched-string BVP numerically. We begin by performing a separation of variables on the wave equation, assuming trial solutions of the form

(3)\[y(x;t) = \tilde{y}(x) \, \exp (-\ii \sigma t),\]

where \(\tilde{y}(x)\) is a function of \(x\) alone. Then, the wave equation reduces to an ordinary differential equation (ODE) for \(\tilde{y}\),

\[\nderiv{\tilde{y}}{x}{2} = - \frac{\sigma^{2}}{c^{2}} \tilde{y}.\]

Discretization

To solve the ODE, we discretize it to establish a set of difference equations. The discretization involves transforming the continuous function \(\tilde{y}(x)\) into a finite set of \(N\) values \(\{\tilde{y}_{1},\tilde{y}_{2},\ldots,\tilde{y}_{N}\}\), representing the function sampled on the discrete spatial grid \(\{x_{1},x_{2},\ldots,x_{N}\}\).

For simplicity let’s assume the grid is uniform, so that

\[x_{j+1} - x_{j} \equiv \Delta x = \frac{L}{N-1} \qquad (1 \leq j \leq N-1).\]

Then, the second derivative of \(\tilde{y}\) can be approximated to second order in \(\Delta x\) as

\[\left. \nderiv{\tilde{y}}{x}{2} \right|_{x=x_{j}} \approx \frac{\tilde{y}_{j+1} - 2 \tilde{y}_{j} + \tilde{y}_{j-1}}{\Delta x^{2}} \qquad (2 \leq j \leq N-1).\]

This allows us to replace the ODE with \(N-2\) difference equations

\[\frac{\tilde{y}_{j+1} - 2 \tilde{y}_{j} + \tilde{y}_{j-1}}{\Delta x^{2}} = - \frac{\sigma^{2}}{c^{2}} \tilde{y}_{j} \qquad (2 \leq j \leq N-1).\]

Together with the two boundary conditions

\[\tilde{y}_{1} = 0 \qquad \tilde{y}_{N} = 0,\]

we thus have a linear system of \(N\) algebraic equations and \(N\) unknowns.

Linear System

To find solutions to the linear system, we first write it in matrix form as

(4)\[\mS \vu = \mathbf{0},\]

where \(\vu\) is the vector with components

\[\begin{split} \vu = \begin{pmatrix} \tilde{y}_{1} \\ \tilde{y}_{2} \\ \vdots \\ \tilde{y}_{N-1} \\ \tilde{y}_{N} \end{pmatrix}\end{split}\]

and the ‘system matrix’ \(\mS\) is an \(N \times N\) tridiagonal matrix with components

\[\begin{split}\mS = \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 1 & \sigma^{2} \tau^{2} - 2 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & \sigma^{2} \tau^{2} - 2 & 1 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 1 \end{pmatrix}.\end{split}\]

Here we’ve introduced

\[\tau \equiv \frac{\Delta x} c\]

as the sound crossing time of a single cell.

Equation (4) is a homogeneous linear system, meaning that it has non-trivial solutions \(\vu\) only when the determinant of \(\mS\) vanishes. With this in mind, we formulate the characteristic equation for the BVP,

\[\Dfunc(\sigma) = 0\]

where \(\Dfunc(\sigma) \equiv \det(\mS)\) is a discriminant function whose roots are the characteristic frequencies (eigenfrequencies) of the stretched-string BVP.

Fig. 1 Plot of the discriminant function \(\Dfunc(\sigma)\) as a function of the frequency \(\sigma\), for the stretched-string BVP with \(N=50\). The orange dots highlight where \(\Dfunc=0\). The function has been scaled so that \(\Dfunc(0) = 1\). (Source)

Fig. 1 plots the discriminant function for the BVP discretized on a spatial grid of \(N=50\) points. The roots (zeros) of the function are highlighted by the orange markers; they fall very close to the values \(\sigma = \pi c/L, 2 \pi c/L, \ldots\) predicted by the analytic solutions.

Scanning for Eigenfrequencies

While Fig. 1 is useful for visualizing \(\Dfunc\), it’s not the best way to find eigenfrequencies. Instead, we can rely on well-established techniques for isolating and refining roots of monovariate functions.

First, we evaluate a set of \(M\) values \(\{\Dfunc_{1},\Dfunc_{2},\ldots,\Dfunc_{M}\}\), representing the discriminant function sampled on the discrete frequency grid \(\{\sigma_{1},\sigma_{2},\ldots,\sigma_{M}\}\). Then, we scan through these data looking for sign changes between adjacent discriminant values. If \(\Dfunc_{i} \Dfunc_{i+1} < 0\), we know that a root of the discriminant function must lie in the interval \((\sigma_{i},\sigma_{i+1})\) — we have bracketed a root. Fig. 2 illustrates the process of bracket scanning for a frequency grid comprising \(M=32\) points, distributed uniformly in \(\sigma\) across the same range as plotted in Fig. 1. This figure highlights five brackets containing the five roots identified previously.

Fig. 2 Plot of the discriminant values \(\{\Dfunc\}\) on the discrete frequency grid \(\{\sigma\}\) (distributed uniformly in \(\sigma\)), for the stretched-string BVP with \(N=50\) and \(M=32\). The orange-haloed segments highlight adjacent points that bracket a root \(\Dfunc=0\). (Source)

Once a bracket is established for a given root, it can be narrowed through a process of iterative refinement until the root is converged upon. There are a variety of well-known root-finding algorithms that perform this refinement; the bisection method is conceptually the simplest, but approaches such as Brent’s method can be much more efficient. For the brackets plotted in Fig. 2, Table 1 compares the eigenfrequencies found using Python’s scipy.optimize.brentq() function, against the analytic values predicted by eqn. (1).

Table 1 Numerical and analytic eigenfrequencies, in units of \(\pi c/L\), for the stretched-string BVP with \(N=50\). (Source)

n	numerical	analytic
1	0.999829	1.000000
2	1.998630	2.000000
3	2.995378	3.000000
4	3.989047	4.000000
5	4.978618	5.000000

Eigenfunction Reconstruction

For each of the eigenfrequencies found, we reconstruct the corresponding eigenfunction by solving the linear system (4). Because \(\det(\mS)\) is now zero, this system is guaranteed to have a non-trivial solution. The solution vector \(\vu\) resides in the null space of \(\mS\), and we can use standard numerical techniques (e.g., singular value decomposition) to evaluate it. Then, the \(j\)’th element of \(\vu\) corresponds to the eigenfunction sampled at the \(j\)’th spatial grid point:

\[(\vu)_{j} = \tilde{y}_{j} \equiv \tilde{y}(x_{j})\]

Fig. 3 Plot of the eigenfunctions \(\tilde{y}\) as a function of spatial coordinate \(x\), for the first three modes of the stretched-string BVP with \(N=50\). The discrete points show the numerical functions, and the solid lines the corresponding analytic functions. In all cases, the eigenfunctions have been normalized to have a maximum \(|\tilde{y}|\) of unity. (Source)

Fig. 3 plots the eigenfunctions found in this way for the first three modes (\(n=1,2,3\)) given in Table 1. Also shown are the corresponding analytic solutions given by eqn. (2). The agreement between the two is good.