# 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

where \(\tilde{y}(x)\) is some 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 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

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

Taken together with the 2 boundary conditions

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

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

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,

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

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:

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.