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

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.

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

 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})$

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.