# The Stretched String Problem

We’ll start our discussion of GYRE by considering the analogous (but much simpler) 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 leads 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(x,t) = B \sin \left( \frac{n \pi x}{L} \right) \exp ( - \ii \sigma t).$

The index $$n$$ uniquely labels the modes, and $$y(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_{k+1} - x_{k} \equiv \Delta x = \frac{L}{N-1} \qquad (1 \leq k \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_{k}} \approx \frac{\tilde{y}_{k+1} - 2 \tilde{y}_{k} + \tilde{y}_{k-1}}{\Delta x^{2}} \qquad (2 \leq k \leq N-1).$

This allows us to replace the ODE with $$N-2$$ difference equations

$\frac{\tilde{y}_{k+1} - 2 \tilde{y}_{k} + \tilde{y}_{k-1}}{\Delta x^{2}} = - \frac{\sigma^{2}}{c^{2}} \tilde{y}_{k} \qquad (2 \leq k \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. 5 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. 5 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_{j} \Dfunc_{j+1} < 0$$, we know that a root of the discriminant function must lie in the interval $$(\sigma_{j},\sigma_{j+1})$$ — we have bracketed a root. Fig. 6 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. 5. 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. 6, Table 1 compares the eigenfrequencies found using Python’s scipy.optimize.brentq() function, against the analytic values predicted by equation (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 $$k$$’th element of $$\vu$$ corresponds to the eigenfunction sampled at the $$k$$’th spatial grid point:

$(\vu)_{k} = \tilde{y}_{k} \equiv \tilde{y}(x_{k})$

Fig. 7 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 equation (2). The agreement between the two is good.