Numerical Solution

Now let’s see how we might go about solving the stretched-string BVP numerically, using an approach very similar to the one implemented in GYRE’s for solving the oscillation equations.

Separation

We begin by performing a separation of variables on the wave equation, assuming trial solutions of the form

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

where \(\tilde{y}(x)\) is some 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 form 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} = \Delta x \equiv \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)\]

Taken together with the 2 boundary conditions

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

we 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

(3)\[\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{L}{N-1} c,\]

representing the sound crossing time of a single cell.

Equation (3) is a homogeneous linear system, meaning that it only has non-trivial solutions \(\vu\) 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.

Plot showing the discriminant function versus frequency

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

Finding Eigenfrequencies

While Fig. 1 is useful for visalizing \(\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 finite 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 inspect the signs of adjacent values \((\Dfunc_{j},\Dfunc_{j+1})\). If these differ, then 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. 2 demonstrates the process of root bracketing for a frequency grid covering the plotted frequency interval with \(M=32\) uniformly spaced points; it highlights five brackets containing the five roots shown previously in Fig. 1.

Plot showing the discriminant function versus frequency, with root brackets indicated

Fig. 2 Plot of the discriminant values \(\{\Dfunc\}\) on the discrete frequency grid \(\{\sigma\}\), for the stretched-string BVP with \(N=50\) and \(M=32\). The orange halos indicate 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 equation (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 find the corresponding eigenfunction by solving the linear system (3). 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})\]
Plot showing eigenfunctions for the first three modes

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

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