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 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 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.
Fig. 2 Plot of the discriminant values \(\{\Dfunc\}\) on the discrete
frequency grid \(\{\sigma\}\) (distributed uniformly in
\(\sigma\)), for the stretched-string BVP with \(N=50\) and
\(M=32\). The orange-haloed segments highlight 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 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 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. In all cases, the eigenfunctions have been
normalized to have a maximum \(|\tilde{y}|\) of
unity. (Source
)
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.