Spectral analysis – the study of the eigenvectors and eigenvalues of linear operators between vector spaces – is one of the most important areas of research in modern applied mathematics, with applications in areas as diverse as mechanics, signal processing, and biology. The spectral analysis of operators involved in partial differential equations is especially important, as it offers keen insight into some of the fundamental laws that make modern engineering possible. Of all these equations, the wave equation in a single dimension is perhaps the simplest and easiest to understand, yet the mathematics that underlie it is both rich and beautiful. In particular, it provides an illuminating platform from which to study spectral theory, for the eigenvalues of the underlying differential operator bear a straightforward physical interpretation as the vibrational frequencies of the string being modeled.

Within spectral theory, there are two broad classes of problems: forward problems and inverse problems. In forward problems, the operator in question is specified, and one is asked to determine its eigenvalues and eigenvectors. In inverse problems, one is instead given information about an operator's eigenstructure and is asked to reconstruct the operator. In the context of the vibrating string problem, the forward problem asks one to find a string's frequencies of vibration given its physical characteristics, while the inverse problem seeks knowledge of the string's physical properties from its frequencies. In the case of the vibrating string, both of these problems have been well studied, and various techniques exist for solving them for strings with both continuous and discrete mass densities.

One possible approach to the spectral analysis of a partial differential equation is to transform it into another form for which the spectral characteristics are known. In the case of the one-dimensional wave equation Equation 1, it is possible to transform it into the Sturm-Liouville equation Equation 6, whose spectral properties have been well-studied (see, e.g., [3]). While such transformations are mathematically very convenient, the original physical interpretation of the problem becomes lost amid changes-of-variables. That is, in order to make physical sense of the solution one must ”untransform" it back into original problem's domain. The focus of our work and of this paper is on this inversion process, specifically for the transformation between the one-dimensional wave equation and its Sturm-Liouville counterpart.

The remainder of this document is divided as follows. In "Problem Setup", we provide a short introduction to the forward and inverse problems under consideration. In "The Transformations", we analyze in detail the transformation used to turn the wave equation into the Sturm-Liouville equation and describe a method for inverting it. In "Numerical Results", we present the results of the main collection of numerical experiments we conducted, which use a particular Sturm-Liouville potential function from [3]. A discussion of the errors in these results is presented in "Error Analysis: Recovering Eigenvalues". "Corresponding q(t) and ρ(x) Functions in Closed Form" contains results for other functions in closed form we studied. Finally, we provide a brief description of the numerical methods we used in "Numerical Methods".

The simplest vibrating string problem is when the string has uniform mass density, ρ(x)=1ρ(x)=1. The wave equation with uniform mass density corresponds to the Sturm-Liouville equation with potential function q(t)=0q(t)=0. For the resulting eigenvalue problem, the eigenvalues are λ={π2,4π2,9π2,16π2,...}λ={π2,4π2,9π2,16π2,...}. What if just one eigenvalue is changed? In chapter 6 of [3] it is shown for μ1<4π2μ1<4π2 the Dirichlet spectra with variable first eigenvalue λ(q)={μ1,4π2,9π2,16π2,...}λ(q)={μ1,4π2,9π2,16π2,...} specify unique, even q(t)q(t)'s given by

where [f,g]=fg'-f'g[f,g]=fg'-f'g, the Wronskian. In [3], a function qq is even if q(1-t)=q(t)q(1-t)=q(t) for t∈[0,1]t∈[0,1].

To illustrate our results in this section, we have set this first eigenvalue, μ1μ1, to 1, 9, and 15. Figure 1 displays the Sturm-Liouville potentials which correspond to these three values of μ1μ1 being substituted into Equation 31.

Figure 1: Sturm-Liouville potentials.

(a) (b) (c)

In solving the initial-value problem described in "Changing From Potential to Density", uniqueness concerns pose the question,"Does a single Sturm-Liouville potential give rise to a single mass density?" As it turns out, no. We see that a single Sturm-Liouville potential can give rise to a set of mass densities, dependent upon the initial conditions ρ(0)ρ(0) and ρ'(0)ρ'(0) and the length of the recovered string, LL. Figure 2 shows a color mapping of recovered string length against various initial conditions. The length of the string cannot be known before the initial-value problem is solved; therefore, in order to obtain these color maps, we fix the initial conditions and obtain the length of the string at the end of the change of variables computation. All string lengths L≥1.25L≥1.25 are indicated on the map as 1.25. We show this result for the transformation of three potential functions, those defined by the aforementioned spectra with μ1=1,9,15μ1=1,9,15. In our case, we desire strings of length one. The curve shown in black on the color maps indicates the set of initial conditions ρ(0)ρ(0) and ρ'(0)ρ'(0) that yield a string of unit length. Alongside this module, in the featured links, we have posted a movie of this color map with μ1μ1 changing from 1 to 15 by 0.1007.

Figure 2: Color map of string lengths recovered from various initial conditions.

(a) (b) (c)

We can determine, from the curve in the color plot, a ρ'(0)ρ'(0) corresponding to any ρ(0)ρ(0) that will yield a string with L=1L=1 for a given μ1μ1. However, we do not use this means. Instead, we choose ρ(0)ρ(0) and use a bisection algorithm, solving the ODE at each iteration, to find ρ'(0)ρ'(0) such that L=1L=1. The bisection algorithm is consistently called with error tolerance τ=10-10τ=10-10. In the results below, we set ρ(0)=1ρ(0)=1 for each μ1μ1-determined potential. Figure 3 shows the mass density functions found in this way, which correspond to the Sturm-Liouville potentials of Figure 1. We have also included a movie in the featured links that shows the mass density with ρ(0)=1ρ(0)=1 and L=1L=1 from potentials determined by μ1μ1 changing from 1 to 37.4 by tenths. (The second eigenvalue in the Dirichlet spectrum for Equation 31 is 4π2≈39.54π2≈39.5. The closer we make μ1μ1 to the second eigenvalue, the harder it is to numerically compute the transformations. Making μ1=37.4μ1=37.4 was as close to 4π24π2 as we were able to compute the transformations.)

Figure 3: ρ(x)ρ(x) when L=1L=1 and ρ(0)=1ρ(0)=1.

(a) (b) (c)

It is also interesting to view a sampling of the isospectral set of mass densities corresponding to a given q(t)q(t). Figures  Figure 4 and  Figure 5 contain four isospectral mass densities, each corresponding to the Sturm-Liouville potential in Equation 31 with μ1μ1 specified as in the other figures and varying values of ρ(0)ρ(0).

Figure 4: Isospectral mass density, μ1=1μ1=1.

(a) (b)

Figure 5: Isospectral mass density functions.

(a) (b)

We have made the claim that the eigenvalues of the Sturm-Liouville operator specified by q(t)q(t) are the same as those of the wave equation specified by the corresponding ρ(x)ρ(x). To demonstrate this, we consider the eigenvalue problem in  Equation 5. The eigenpairs can be numerically obtained via spectral methods. We build a differential matrix, DD, over a Chebyshev discretization of the interval and in MATLAB solve the linear system:

A short description of spectral discretization is in "Numerical Methods".

For the specific qq in Equation 31, we know what the eigenvalues are supposed to be, so we can calculate them and analyze the error. Before doing so, it is interesting to consider the eigenfunctions for several ρ(x)ρ(x) obtained from different μ1μ1 values. Figure 6 shows the first two eigenfunctions obtained from the numerically calculated mass density functions. When μ1=9μ1=9, we see that the first two eigenfunctions resemble φ1(x)=sin(πx)φ1(x)=sin(πx) and φ2(x)=sin(2πx)φ2(x)=sin(2πx). This is to be expected because μ1=9μ1=9 is fairly close to π2≈9.9π2≈9.9. (When μ1=π2μ1=π2, we obtain q(t)=0q(t)=0 from  Equation 31 and therefore a string of uniform mass density, which has as its first two eigenfunctions sin(πx)sin(πx) and sin(2πx)sin(2πx).) Equivalently, potentials with μ1μ1 farthest from π2π2 show eigenfunction behavior most deviant from sin(πx)sin(πx) and sin(2πx)sin(2πx). Interested readers can view a movie of the eigenfunctions changing as μ1μ1 changes from 1 to 37.4 by tenths. The behavior of the eigenfunctions as μ1μ1 increases changes drastically when μ1≥35μ1≥35 due to the first and second eigenvalues becoming close. In the following section, we discuss the error in numerically calculating the eigenvalues.

Figure 6: The first two eigenfunctions recovered from the mass densities shown in Figure 3. The first eigenfunction is shown in blue, the second in red.

The sources of error in our numerical scheme arise in part because of the highly constrained nature of the ODEs. Unless otherwise specified, the integration was completed with step size 10-410-4. In addition to this error, we are required to interpolate between nodes to obtain the values of ρ(x)ρ(x) at the points in the Chebyshev grid. (Unless otherwise specified, all computations were executed on a Chebyshev grid of size N=256N=256.) However, it should be noted that in all runs, we have a minimum of 104104 nodes between 0 and 1. With a grid so fine, even linear interpolation to the Chebyshev points yields results with tolerable error. Because, in our case, the potential function q(t)q(t) was originally specified by its Dirichlet spectrum, we know the exact eigenvalues that the numerical solver should return. Figure 7 shows the first five expected eigenvalues alongside those we recovered numerically. The error is visibly small, but indistinguishable from case to case. Figure 8 shows the error in the first eigenvalue recovered from transformations of potential functions specified by various μ1μ1. The error is smallest near μ1=π2μ1=π2, which is to be expected because this is the simplest case. As μ1μ1 gets further away from π2π2, the error increases. We also see that decreasing the time step in the ODE solver by a single order of magnitude reduces the error by a single order of magnitude. In this way, we have verified first-order convergence, which is consistent with the expected results of the forward Euler scheme used.

Figure 7: The first five eigenvalues of q(t)q(t) (blue circles) and those numerically recovered from the transformed equivalent, ρ(x),x∈[0,1]ρ(x),x∈[0,1] (green 'x').

Figure 8: The error in the first eigenvalue recovered from the transformation of q(t)q(t) with λ(q)={μ1,4π2,9π2,16π2,...}λ(q)={μ1,4π2,9π2,16π2,...} to ρ(x),x∈[0,1]ρ(x),x∈[0,1] with varying step size.

The following tables display the error in the first 10 recovered eigenvalues. Listed are the eigenvalues numerically recovered from q(t)q(t) and those recovered from its transformed counterpart, ρ(x)ρ(x). As expected, those recovered from q(t)q(t) are accurate to over ten decimal digits. This is a consequence of the accuracy of the spectral method with sufficient Chebyshev points. The error in the eigenvalues recovered from ρ(x)ρ(x) reflects the numerical error in our ODE solver. However, the error is tolerable. In addition, we see that the error in the transformation of the q(t)q(t) defined by μ1=9μ1=9 is less than for the other choices of μ1μ1 listed in the tables.

For the majority of our research, the focus was on the potential function Equation 31 from chapter 6 of [3]. In the next section, some results will be discussed for other q(t)q(t) or ρ(x)ρ(x) in closed form.

The potential function q(t)q(t) given by Equation 31 from chapter 6 of [3], after taking the second derivative of the Wronskian, becomes

where f(t)= sin μ1(1-t)- sin μ1t sin μ1f(t)= sin μ1(1-t)- sin μ1t sin μ1 and f'(t)=-μ1 cos μ1(1-t)+ cos μ1t sin μ1f'(t)=-μ1 cos μ1(1-t)+ cos μ1t sin μ1. To compare, in the rest of this section, some results will be discussed for other q(t)q(t) or ρ(x)ρ(x) which can be written much more simply.