# 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

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

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,

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

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

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

and so

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

and the corresponding eigenfunctions are

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

where \(\tilde{y}(x)\) is a function of \(x\) alone. Then, the wave equation reduces to an ordinary differential equation (ODE) for \(\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

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

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

Together with the two boundary conditions

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

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

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

Here we’ve introduced

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,

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

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.

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).

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:

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.