# Connexions

You are here: Home » Content » The Art of the PFUG » There and Back Again - An Exploration of the Liouville Transformation and Its Inverse

### Lenses

What is a lens?

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

#### Affiliated with (What does "Affiliated with" mean?)

This content is either by members of the organizations listed or about topics related to the organizations listed. Click each link to see a list of all content affiliated with the organization.
• Rice Digital Scholarship

This collection is included in aLens by: Digital Scholarship at Rice University

Click the "Rice Digital Scholarship" link to see all content affiliated with them.

#### Also in these lenses

• Lens for Engineering

This collection is included inLens: Lens for Engineering
By: Sidney Burrus

Click the "Lens for Engineering" link to see all content selected in this lens.

### Recently Viewed

This feature requires Javascript to be enabled.

Inside Collection (Book):

Book by: Steven J. Cox. E-mail the author

# There and Back Again - An Exploration of the Liouville Transformation and Its Inverse

Module by: John Vogelgesang, Sharmaine Jennings, Anthony Austin. E-mail the authorsEdited By: John Vogelgesang, Sharmaine Jennings, Anthony Austin

Summary: This report summarizes work done as part of the Physics of Strings PFUG under Rice University's VIGRE program. VIGRE is a program of Vertically Integrated Grants for Research and Education in the Mathematical Sciences under the direction of the National Science Foundation. A PFUG is a group of Postdocs, Faculty, Undergraduates and Graduate students formed around the study of a common problem. This module describes a method for recovering the distribution of mass along a one-dimensional non-uniform string from the string's eigenvalues. The eigenvalues are used to construct the string's Sturm-Liouville potential function, and then the mass density is obtained from this potential by means of inverting the Liouville transformation.

## Introduction

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

## Problem Setup

### The Forward Problem

We begin by modeling the behavior of a string of length LL with spatially-varying mass density ρ(x)>0ρ(x)>0 for x[0,L]x[0,L]. For a long string experiencing small vibration, let u(x,t)u(x,t) be the vertical displacement of the string at time t0t0. This u(x,t)u(x,t) satisfies the 1-D wave equation

ρ ( x ) u t t ( x , t ) = u x x ( x , t ) , 0 x L , ρ ( x ) u t t ( x , t ) = u x x ( x , t ) , 0 x L ,
(1)

subject to the Dirichlet boundary conditions

u ( 0 , t ) = u ( L , t ) = 0 , t 0 u ( 0 , t ) = u ( L , t ) = 0 , t 0
(2)

where utt(x,t)=2ut2utt(x,t)=2ut2. Assuming the variables can be separated, we write u(x,t)=ψ(t)y(x)u(x,t)=ψ(t)y(x) and after making the appropriate substitutions, obtain

ρ ( x ) ψ ' ' ( t ) y ( x ) = ψ ( t ) y ' ' ( x ) , ρ ( x ) ψ ' ' ( t ) y ( x ) = ψ ( t ) y ' ' ( x ) ,
(3)

hence

1 ρ ( x ) y ' ' ( x ) y ( x ) = ψ ' ' ( t ) ψ ( t ) . 1 ρ ( x ) y ' ' ( x ) y ( x ) = ψ ' ' ( t ) ψ ( t ) .
(4)

As the left side depends only on xx and the right side depends only on tt, we conclude that both sides must be equal to the same constant, which we denote by -λ-λ. Equating the left side with -λ-λ and rearranging yields

- y ' ' ( x ) = λ ρ ( x ) y ( x ) - y ' ' ( x ) = λ ρ ( x ) y ( x )
(5)

with boundary conditions y(0)=y(L)=0y(0)=y(L)=0. The constant λλ is said to be an eigenvalue of the boundary-value problem Equation 1-Equation 2, and y(x)y(x) is its corresponding eigenfunction. Finding all λλ such that Equation 5 has a solution satisfying the boundary conditions is called the forward eigenvalue problem, and the collection of all such λλ is called the problem's spectrum. Physically, the forward problem asks one to determine a string's natural frequencies from a description of its physical composition.

The inverse problem will do the opposite, asking for information about the string from its eigenvalues.

### The Inverse Problem

In contrast to the forward problem, the inverse eigenvalue problem asks if one can recover information about a string's physical composition from knowledge of its frequencies. Mathematically, this amounts to asking if we can determine the mass density function ρ(x)ρ(x) given the spectrum of the problem Equation 1-Equation 2.

In determining physical properties of objects, the inverse eigenvalue problem is involved in answering various questions. For a vibrating string one question is, "What is the isospectral set, the mass density of all strings of a given length with a given set of frequencies?" The inverse eigenvalue problem could be used in designing a musical instrument that will produce a desired sound. Also, when a particular object is being designed and it is known that the object will be experiencing vibrations of certain frequencies, the inverse problem can be used to ensure the object will not resonate at these frequencies. Additionally, the inverse problem is involved in structural composition analysis. By collecting acoustic data for a structure, a bridge or building, for example, it can be determined if there are "cracks" or "faults" in the object's components. This allows for analysis of a structure without dismantling the structure.

The inverse problem can be attacked by many means, including several numerical schemes. In order to apply the techniques described in [3], it is necessary to transform Equation 5 into the Sturm-Liouville equation

- w ' ' + q w = λ w , - w ' ' + q w = λ w ,
(6)

where qq is called the potential. It is the transformation between these equations that captures the focus of our research and of this paper.

## The Transformations

### Changing From Mass Density to Sturm-Liouville Potential

To transform Equation 5 into the Sturm-Liouville equation Equation 6, it is necessary to employ a change of variables that is due to Joseph Liouville [2]. Liouville's transformation requires ρC2ρC2[0, L]. If that condition is satisfied, the transformation involves a change of dependent and independent variables. The independent variable xx is transformed into a new independent variable tt according to

t = 0 x ρ ( α ) d α . t = 0 x ρ ( α ) d α .
(7)

(NB: This tt is completely unrelated to the tt variable used to denote the dependence on time in Equation 1.) The new dependent variable w(t)w(t) is given by

w ( t ) = y ( x ) ρ 1 / 4 ( x ) , w ( t ) = y ( x ) ρ 1 / 4 ( x ) ,
(8)

where the variable xx on the right-hand side of this equation is meant to be interpreted as a function of tt whose definition is given implicitly by the integral Equation 7.

Taking the derivative once,

d w d t d t d x = 1 4 y ( x ) [ ρ ( x ) ] - 3 / 4 d ρ d x + d y d x [ ρ ( x ) ] 1 / 4 d w d t d t d x = 1 4 y ( x ) [ ρ ( x ) ] - 3 / 4 d ρ d x + d y d x [ ρ ( x ) ] 1 / 4
(9)

(suppressing the argument of x for now)

d w d t ρ 1 / 2 = 1 4 y ρ - 3 / 4 ρ ' + y ' ρ 1 / 4 d w d t ρ 1 / 2 = 1 4 y ρ - 3 / 4 ρ ' + y ' ρ 1 / 4
(10)

Taking the derivative again,

d 2 w d t 2 d t d x ρ 1 / 2 + 1 2 d w d t ρ - 1 / 2 ρ ' = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' + y ' ' ρ 1 / 4 d 2 w d t 2 d t d x ρ 1 / 2 + 1 2 d w d t ρ - 1 / 2 ρ ' = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' + y ' ' ρ 1 / 4
(11)

Substitution and simplifying terms leads to

d 2 w d t 2 ρ + 1 2 d w d t ρ - 1 / 2 ρ ' = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' + y ' ' ρ 1 / 4 d 2 w d t 2 ρ + 1 2 d w d t ρ - 1 / 2 ρ ' = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' + y ' ' ρ 1 / 4
(12)
d 2 w d t 2 ρ + ρ ' 2 ρ 1 / 2 y ρ ' 4 ρ 5 / 4 + y ' ρ 1 / 4 = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' + y ' ' ρ 1 / 4 d 2 w d t 2 ρ + ρ ' 2 ρ 1 / 2 y ρ ' 4 ρ 5 / 4 + y ' ρ 1 / 4 = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' + y ' ' ρ 1 / 4
(13)
d 2 w d t 2 ρ + ρ ' 2 ρ 1 / 2 y ρ ' 4 ρ 5 / 4 + y ' ρ 1 / 4 = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' - λ y ρ 5 / 4 d 2 w d t 2 ρ + ρ ' 2 ρ 1 / 2 y ρ ' 4 ρ 5 / 4 + y ' ρ 1 / 4 = 1 2 y ' ρ - 3 / 4 ρ ' - 3 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' - λ y ρ 5 / 4
(14)
d 2 w d t 2 ρ = - 5 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' - λ y ρ 5 / 4 d 2 w d t 2 ρ = - 5 16 y ρ - 7 / 4 ( ρ ' ) 2 + 1 4 y ρ - 3 / 4 ρ ' ' - λ y ρ 5 / 4
(15)
d 2 w d t 2 = - 5 16 y ρ - 11 / 4 ( ρ ' ) 2 + 1 4 y ρ - 7 / 4 ρ ' ' - λ y ρ 1 / 4 d 2 w d t 2 = - 5 16 y ρ - 11 / 4 ( ρ ' ) 2 + 1 4 y ρ - 7 / 4 ρ ' ' - λ y ρ 1 / 4
(16)
d 2 w d t 2 = - 5 16 y ρ 1 / 4 ρ - 3 ( ρ ' ) 2 + 1 4 y ρ 1 / 4 ρ - 2 ρ ' ' - λ y ρ 1 / 4 d 2 w d t 2 = - 5 16 y ρ 1 / 4 ρ - 3 ( ρ ' ) 2 + 1 4 y ρ 1 / 4 ρ - 2 ρ ' ' - λ y ρ 1 / 4
(17)
d 2 w d t 2 = - 5 16 w ρ - 3 ( ρ ' ) 2 + 1 4 w ρ - 2 ρ ' ' - λ w d 2 w d t 2 = - 5 16 w ρ - 3 ( ρ ' ) 2 + 1 4 w ρ - 2 ρ ' ' - λ w
(18)
- d 2 w d t 2 + 1 4 ρ - 2 ρ ' ' - 5 16 ρ - 3 ( ρ ' ) 2 w = λ w - d 2 w d t 2 + 1 4 ρ - 2 ρ ' ' - 5 16 ρ - 3 ( ρ ' ) 2 w = λ w
(19)

This equation is now in the form of the Sturm-Liouville equation

- d 2 w d t 2 + q ( t ) w = λ w , - d 2 w d t 2 + q ( t ) w = λ w ,
(20)

where the Sturm-Liouville Dirichlet potential function, q(t)q(t), is given by

q ( t ) = ρ ' ' ( x ) 4 ρ 2 ( x ) - 5 ρ ' ( x ) 2 16 ρ 3 ( x ) . q ( t ) = ρ ' ' ( x ) 4 ρ 2 ( x ) - 5 ρ ' ( x ) 2 16 ρ 3 ( x ) .
(21)

As in Equation 8 above, the variable xx on the right side of this expression is to be interpreted as a function of tt defined by Equation 7.

This section describes what is needed to change from the wave equation with a mass density, ρρ, to a Sturm-Liouville equation with a potential, qq. The next section describes the reverse, changing from the Sturm-Liouville equation with a potential, qq, to a wave equation with a mass density, ρρ.

### Changing From Potential to Density

Converting a given Dirichlet potential function q(t)q(t) into corresponding mass density ρ(x)ρ(x) amounts to solving the nonlinear ordinary differential equation, or ODE, in Equation 21 for ρρ; however, one must be careful. Using the integral Equation 7 and abusing notation a bit, it is tempting to write t=t(x)t=t(x) and arrive at the ODE

q t ( x ) = ρ ' ' ( x ) 4 ρ 2 ( x ) - 5 ρ ' ( x ) 2 16 ρ 3 ( x ) , q t ( x ) = ρ ' ' ( x ) 4 ρ 2 ( x ) - 5 ρ ' ( x ) 2 16 ρ 3 ( x ) ,
(22)

which one would then attempt to solve numerically by selecting some initial conditions for ρ(x)ρ(x), picking a grid in the xx-domain and applying one's ODE solver of choice, using some quadrature scheme to estimate t(x)t(x) by approximating Equation 7. This approach, however, suffers from a somewhat subtle flaw: it presumes that one already knows the domain of definition of ρ(x)ρ(x), i.e., that the length LL of the string can be known a priori.

Observe that the ODE makes sense only over the domain of definition of the tt variable, while the function ρρ makes sense only over the domain of definition of the xx variable. Our goal was to learn about the qq in Chapter 6 of [3], shown in Equation 31. Since the function qq is known, we have knowledge of the domain of tt. In principle, we can obtain the corresponding domain of xx by inverting the integral Equation 7, but this requires knowledge of ρρ as a function of xx, which is what we are trying to find in the first place. We need a way to obtain this domain information from other quantities we can compute.

Our answer to this problem is the following. Treating Equation 7 as an equation that defines xx implicitly as a function of tt, we differentiate both sides with respect to tt to obtain

1 = ρ 1 / 2 x ( t ) x ' ( t ) , 1 = ρ 1 / 2 x ( t ) x ' ( t ) ,
(23)

which can be rearranged to give

x ' ( t ) = ρ - 1 / 2 x ( t ) . x ' ( t ) = ρ - 1 / 2 x ( t ) .
(24)

Thus, we obtain a differential equation for x(t)x(t). Furthermore, we have an initial condition: x(0)=0x(0)=0, since x=0x=0 and t=0t=0 correspond to each other by Equation 7. We can therefore obtain the xx-values that correspond to given tt-values, provided that we can compute ρ-1/2x(t)ρ-1/2x(t) for arbitrary tt.

We accomplish this by changing the equation Equation 21 into an ODE for σ(t)=(ρx)(t)=ρx(t)σ(t)=(ρx)(t)=ρx(t). By the chain rule,

σ ' ( t ) = ρ ' x ( t ) · x ' ( t ) . σ ' ( t ) = ρ ' x ( t ) · x ' ( t ) .
(25)

Since x'(t)=σ-1/2(t)x'(t)=σ-1/2(t), we obtain

σ ' ( t ) = ρ ' x ( t ) σ ( t ) . σ ' ( t ) = ρ ' x ( t ) σ ( t ) .
(26)

Applying the quotient rule produces

σ ' ' ( t ) = σ ( t ) ρ ' ' x ( t ) x ' ( t ) - 1 2 ρ ' x ( t ) σ - 1 / 2 ( t ) σ ' ( t ) σ ( t ) = σ ( t ) ρ ' ' x ( t ) σ - 1 / 2 ( t ) - 1 2 σ ' ( t ) 2 σ ( t ) = ρ ' ' x ( t ) σ ( t ) - σ ' ( t ) 2 2 σ ( t ) . σ ' ' ( t ) = σ ( t ) ρ ' ' x ( t ) x ' ( t ) - 1 2 ρ ' x ( t ) σ - 1 / 2 ( t ) σ ' ( t ) σ ( t ) = σ ( t ) ρ ' ' x ( t ) σ - 1 / 2 ( t ) - 1 2 σ ' ( t ) 2 σ ( t ) = ρ ' ' x ( t ) σ ( t ) - σ ' ( t ) 2 2 σ ( t ) .
(27)

Hence,

ρ ' ' x ( t ) = σ ( t ) σ ' ' ( t ) + 1 2 σ ' ( t ) 2 . ρ ' ' x ( t ) = σ ( t ) σ ' ' ( t ) + 1 2 σ ' ( t ) 2 .
(28)

We can substitute these expressions for ρ'x(t)ρ'x(t) and ρ''x(t)ρ''x(t) in equation Equation 21 to obtain

q ( t ) = σ ( t ) σ ' ' ( t ) + 1 2 σ ' ( t ) 2 4 σ ( t ) 2 - 5 σ ' ( t ) 2 σ ( t ) 16 σ ( t ) 3 , q ( t ) = σ ( t ) σ ' ' ( t ) + 1 2 σ ' ( t ) 2 4 σ ( t ) 2 - 5 σ ' ( t ) 2 σ ( t ) 16 σ ( t ) 3 ,
(29)

which may be simplified to

q ( t ) = σ ' ' ( t ) 4 σ ( t ) - 3 16 σ ' ( t ) σ ( t ) 2 . q ( t ) = σ ' ' ( t ) 4 σ ( t ) - 3 16 σ ' ( t ) σ ( t ) 2 .
(30)

This is a differential equation that we may solve for σ(t)σ(t) over a grid that coincides with the domain of qq. We can then use the ODE for x'(t)x'(t), Equation 24, to obtain the values of xx at which ρρ takes on the values of σ(t)σ(t).

We have shown that it is possible to convert a mass density function, ρ(x)ρ(x), x[0,L]x[0,L], in the wave equation to a potential function, q(t)q(t), t[0,1]t[0,1], in the Sturm-Liouville equation and discussed our method to numerically convert from qq to ρρ. In the next section, the results from numerically reverting this change of variables, arriving at a mass density function for a string of unit length from a specific Sturm-Liouville potential, are described.

## Numerical Results

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

q ( t ) = - 2 d 2 d t 2 log sin μ 1 ( 1 - t ) - sin μ 1 t sin μ 1 , sin π t q ( t ) = - 2 d 2 d t 2 log sin μ 1 ( 1 - t ) - sin μ 1 t sin μ 1 , sin π t
(31)

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.

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

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π239.54π239.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.)

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

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:

- D 2 y ( x ) = λ ρ ( x ) y ( x ) . - D 2 y ( x ) = λ ρ ( x ) y ( x ) .
(32)

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 π29.9π29.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 μ135μ135 due to the first and second eigenvalues becoming close. In the following section, we discuss the error in numerically calculating the eigenvalues.

## Error Analysis: Recovering Eigenvalues

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.

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.

 μ 1 = 1 μ 1 = 1 True λλ Recovered λ(q(t))λ(q(t)) Error from q(t)q(t) Recovered λ(ρ(x))λ(ρ(x)) Error from ρ(x)ρ(x) 1.0000000 1.0000000 2.5930147e-011 1.0054709 5.4708902e-003 39.4784176 39.4784176 1.0608403e-011 39.5200555 4.1637860e-002 88.8264396 88.8264396 1.5930368e-011 88.9392066 1.1276694e-001 157.9136704 157.9136704 1.4694024e-011 157.9977790 8.4108586e-002 246.7401100 246.7401100 2.2168933e-012 247.2588167 5.1870663e-001 355.3057584 355.3057584 1.6484591e-011 354.0659603 1.2397981e+000 483.6106157 483.6106157 8.5265128e-013 486.5704124 2.9597967e+000 631.6546817 631.6546817 5.1159077e-012 620.2338813 1.1420800e+001 799.4379565 799.4379565 9.7770680e-012 816.7159848 1.7278028e+001 986.9604401 986.9604401 2.6602720e-011 973.9089831 1.3051457e+001
 μ 1 = 9 μ 1 = 9 True λλ Recovered λ(q(t))λ(q(t)) Error from q(t)q(t) Recovered λ(ρ(x))λ(ρ(x)) Error from ρ(x)ρ(x) 9.0000000 9.0000000 2.3881341e-011 9.0000369 3.6905319e-005 39.4784176 39.4784176 2.0996538e-011 39.4791540 7.3639246e-004 88.8264396 88.8264396 2.1927349e-011 88.8280954 1.6558273e-003 157.9136704 157.9136704 1.7990942e-011 157.9166007 2.9303014e-003 246.7401100 246.7401100 4.2064130e-012 246.7446791 4.5690754e-003 355.3057584 355.3057584 2.1827873e-011 355.3123306 6.5721222e-003 483.6106157 483.6106157 2.1032065e-011 483.6195551 8.9394021e-003 631.6546817 631.6546817 3.3082870e-011 631.6663524 1.1670772e-002 799.4379565 799.4379565 1.3528734e-011 799.4527229 1.4766388e-002 986.9604401 986.9604401 3.4106051e-013 986.9786662 1.8226098e-002
 μ 1 = 15 μ 1 = 15 True λλ Recovered λ(q(t))λ(q(t)) Error from q(t)q(t) Recovered λ(ρ(x))λ(ρ(x)) Error from ρ(x)ρ(x) 15.0000000 15.0000000 8.2795992e-012 14.9930764 6.9236370e-003 39.4784176 39.4784176 1.5845103e-012 39.4741052 4.3124162e-003 88.8264396 88.8264396 6.0538241e-012 88.8165421 9.8974784e-003 157.9136704 157.9136704 8.0433438e-012 157.8982531 1.5417339e-002 246.7401100 246.7401100 2.8421709e-014 246.7173201 2.2789901e-002 355.3057584 355.3057584 1.1766588e-011 355.2739142 3.1844196e-002 483.6106157 483.6106157 2.8762770e-011 483.5680542 4.2561441e-002 631.6546817 631.6546817 8.7538865e-012 631.5997459 5.4935770e-002 799.4379565 799.4379565 3.4106051e-011 799.3689927 6.8963814e-002 986.9604401 986.9604401 2.3192115e-011 986.8757937 8.4646413e-002

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.

## Corresponding q(t)q(t) and ρ(x)ρ(x) Functions 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

q ( t ) = - 2 - π 3 cos ( π t ) f ( t ) - π 2 sin ( π t ) f ' ( t ) + π μ 1 cos ( π t ) f ( t ) + μ 1 sin ( π t ) f ' ( t ) π cos ( π t ) f ( t ) - sin ( π t ) f ' ( t ) = - - π 2 sin ( π t ) f ( t ) + μ 1 sin ( π t ) f ( t ) 2 - π cos ( π t ) f ( t ) - sin ( π t ) f ' ( t ) 2 q ( t ) = - 2 - π 3 cos ( π t ) f ( t ) - π 2 sin ( π t ) f ' ( t ) + π μ 1 cos ( π t ) f ( t ) + μ 1 sin ( π t ) f ' ( t ) π cos ( π t ) f ( t ) - sin ( π t ) f ' ( t ) = - - π 2 sin ( π t ) f ( t ) + μ 1 sin ( π t ) f ( t ) 2 - π cos ( π t ) f ( t ) - sin ( π t ) f ' ( t ) 2
(33)

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.

### qq and ρρ that satisfy the Liouville Transformation

The potential function and the corresponding mass density function

q ( t ) = - 3 ( 1 + 4 t ) 2 and ρ ( x ) = ( 2 x + 1 ) 2 , q ( t ) = - 3 ( 1 + 4 t ) 2 and ρ ( x ) = ( 2 x + 1 ) 2 ,
(34)

respectively, satisfy Equation 21, the Liouville transformation. The graphs of q(t)q(t) and ρ(x)ρ(x) are given in Figure 9 and Figure 10, respectively. In comparing these graphs with the results from the potential function Equation 31, one difference to note is that the qq in Equation 34 has a domain of t[0,2]t[0,2]. This potential is not even, over any interval of positive tt values, which is another difference from the potential from [3]. As a result of these differences, it is not surprising that the ρρ in Equation 34 has a different graph from all the computed ρρ from Equation 31 for the various values of μ1μ1, but there are some similarities. This ρρ in closed form, has ρ(0)=1ρ(0)=1 as we required for our results in "Numerical Results". The graph of this ρρ in Figure 10 has another similarity to the computed ρρ when μ1=17.1μ1=17.1 because that graph also shows ρ(1)=9ρ(1)=9, however, the shape of the graphs between those two endpoints is quite different.

The spacing of the eigenvalues produced by Equation 34, calculated by spectral discretization, can be seen in Figure 11. For comparison, in addition to the computed eigenvalues from these equations, the eigenvalues for the string of uniform density are graphed as small black dots in Figure 12. The expected eigenvalues from Equation 31 are the same as the eigenvalues for the string of uniform density, except for the first one. The eigenvalues from Equation 34 show the same quadratic behavior as is seen in Figure 7, but these eigenvalues increase very slowly compared to the eigenvalues for the uniform mass density.

Next is an example of a density function in closed form from Euler.

### Example from Euler

Leonhard Euler found the wave equation has a general solution of

y = 1 + x α Φ x x + x α + c 0 t + Ψ x 1 + x α - c 0 t y = 1 + x α Φ x x + x α + c 0 t + Ψ x 1 + x α - c 0 t
(35)

where c0Tσ0c0Tσ0 when the "line density" is given by σ(x)=σ01+xα4σ(x)=σ01+xα4[5]. Euler's analysis showed the frequencies of this string would be vk=k2L1+LαTσ0vk=k2L1+LαTσ0. In order to produce figures which would compare easily to the results from the previous section, the following constants were chosen: σ0=3σ0=3 and α=12α=12. The resultant graph of σ(x)σ(x), which is very different from the other mass density graphs in this paper, can be seen in Figure 13, and the first 256 frequencies are graphed in Figure 14.

The relationship between eigenvalues, λkλk, and frequencies, vkvk, is given by λk=2πvk2λk=2πvk2. Using this relationship, in Figure 15, the exact eigenvalues for k=1k=1 to 256 have been graphed as blue circles. To compare the results from this set of equations to the others in this module, the eigenvalues were computed for this σ(x)σ(x) using spectral discretization with N=256N=256. These computed eigenvalues are also graphed in Figure 15, each marked with a green 'x'. The error between the eigenvalues from the frequency function and the computed eigenvalues is graphed in Figure 16, which shows the first 100 eigenvalues are computed with reasonable accuracy. For comparison, in addition to the exact and computed eigenvalues from Euler's equations, the eigenvalues for the string of uniform density are graphed as black dots in Figure 17. The eigenvalues from Euler's equations still show the same quadratic behavior as the other examples we have seen, but these eigenvalues increase very rapidly compared to the eigenvalues for the uniform mass density.

In computing the results which are reported starting in "Numerical Results", various techniques were used. Some techniques were successful and other techniques failed to give the desired or expected results, which led to alternate techniques being tried. A short description of our methods, including a description of spectral discretization, is given in the next section.

## Numerical Methods

When we began gathering the data included here in our results, we started by programming the q(t)q(t) in Equation 31 from [3] into MATLAB. Our next goal was to find the ρ(x)ρ(x) from these potential functions. We started by assuming x[0,1]x[0,1] and used a uniform grid to solve the ODE from the Liouville transformation Equation 21 with forward Euler and the change of variables for tt Equation 7 with composite trapezoid quadrature. Briefly, forward Euler gives a numerical approximation to the solution of the initial value problem

a ' = f ( a , b ) , a ( b 0 ) = a 0 a ' = f ( a , b ) , a ( b 0 ) = a 0
(36)

through the formula

a n + 1 = a n + h f ( a n , b n ) a n + 1 = a n + h f ( a n , b n )
(37)

where h=bn+1-bnh=bn+1-bn[1]. In order to facilitate finding the eigenvalues, we attempted to incorporate a Chebyshev grid for xx into the uniform grid. This is when we became aware something was not working properly. In an attempt to find the problem, we output the value of tt at each step, and we realized for some choices of μ1μ1 and initial conditions, when xx reached 1, tt might be much less than 1 or even greater than 1. We had two choices at this point, we could choose values of μ1μ1 and initial conditions and find the necessary value of xx which would result in tt being within some tolerance of 1 or we could solve the ODE in tt instead of xx. This second method is how we chose to continue our exploration.

With a uniform grid over t[0,1]t[0,1] and given μ1μ1, ρ(0)ρ(0), and ρ'(0)ρ'(0) we used forward Euler to solve the ODE Equation 30 and then implemented the change of variables. This resulted in strings of various lengths depending on the choice of μ1μ1, ρ(0)ρ(0), and ρ'(0)ρ'(0). These are illustrated in Figure 2.

We made the choice to restrict to strings of unit length and ρ(0)=1ρ(0)=1. To facilitate this choice, we included a bisection solver to find the necessary value of ρ'(0)ρ'(0) for a given μ1μ1 that results in a string of length 1. After this determination is made, the resulting graphs are illustrated in Figure 3.

For a given μ1μ1 by changing the choice of ρ(0)ρ(0), with the necessary value of ρ'(0)ρ'(0) being determined to result in a unit string, the isospectral graphs, graphs with the same Dirichlet spectrum, are produced (seen in Figure 4).

To compute the eigenvalues and eigenfunctions, we wanted to use spectral discretization, but because our method solves the ODE in tt, we could not choose at which xx values we would have values for ρ(x)ρ(x). To overcome this, we used the MATLAB function "interp1" and used a cubic spline interpolation to find ρρ at the Chebyshev grid x values from the values we did have from solving the ODE in tt. Then, after adjusting the differential operator matrix, DD, and ρρ matrix, BB, for the Dirichlet boundary conditions, we used the MATLAB function "eig" to solve for MM, the diagonal matrix of generalized eigenvalues, and VV, a full matrix of corresponding eigenvectors, in the generalized eigenvalue problem, D2*V=M*B*VD2*V=M*B*V (from -y''(x)=λρ(x)y(x)-y''(x)=λρ(x)y(x)). The eigenvalues were sorted, and the eigenvectors were sorted correspondingly, to produce the figures in this paper.

In an effort to improve accuracy, Heun's method was implemented instead of forward Euler. In Heun's method, the formula for the numerical approximation to the solution of Equation 36 is

a n + 1 = a n + h 2 f ( a n , b n ) + f a n + h f ( a n , b n ) , b n + 1 . a n + 1 = a n + h 2 f ( a n , b n ) + f a n + h f ( a n , b n ) , b n + 1 .
(38)

where h=bn+1-bnh=bn+1-bn[1]. Additionally, RK4, the fourth order Runge-Kutta method, was implemented through "ode45" in MATLAB. However, these methods did not produce results that were much better than the simpler method of forward Euler (and in some cases, took much more time), so we left forward Euler in our code and reported those results in this paper.

### Spectral Discretization

To numerically compute the eigenvalues for the various operators resulting from the various qq and ρρ functions discussed in this paper, there are several techniques. One technique to approximate a second derivative is to use a uniform grid, 0=x0<x1<...<xn=10=x0<x1<...<xn=1 with xj+1-xj=hxj+1-xj=h for j=0,1,...,n-1j=0,1,...,n-1, and a finite difference which leads to

u ' ' ( x j ) u ( x j + 1 ) - 2 u ( x j ) + u ( x j - 1 ) h 2 u ' ' ( x j ) u ( x j + 1 ) - 2 u ( x j ) + u ( x j - 1 ) h 2
(39)

for j=1,...,n-1j=1,...,n-1 which is O(h2)O(h2) accurate. Another technique is spectral discretization. In this technique, instead of a uniform grid, Chebyshev points are used, which are given by

x j = 1 2 1 - cos j π n x j = 1 2 1 - cos j π n
(40)

for j=0,1,...,nj=0,1,...,n. These points are clustered near the endpoints of the interval. These Chebyshev points are used to form a polynomial, pn(x)pn(x), which is an interpolant of u(x)u(x) so that uu(xj)pn''(xj)(xj)pn''(xj). The MATLAB program we used to create the matrix that approximates the differential operator and the Chebyshev grid for xx is based on "cheb.m" in [4].

The advantage of spectral discretization over the finite difference technique is that with the same number of points in the grid, spectral discretization produces results that are much more accurate. Or, it is possible to achieve the same accuracy as finite difference with fewer points using spectral discretization. In Figure 15, spectral discretization was used with N=256N=256. Even though there are infinitely many eigenvalues, this technique's results only approximate the first NN eigenvalues. Since for this problem, it is known exactly what the eigenvalues should be, the error was plotted in Figure 16, and about the first 100 eigenvalues are accurate to 10-710-7 or better. To compare, using finite difference with N=256N=256 and a uniform grid, h=1/256h=1/256 and O(h2)O(h2) accuracy means h21.5×10-5h21.5×10-5. Spectral discretization produces more digits of accuracy than finite difference with the same number of grid points.

## Summary

In this paper, we started by describing Liouville's transformation to change the mass density of the wave equation to the potential of the Sturm-Liouville equation [2]. However, our intention was to study a specific qq from [3], so it was necessary to reverse that transformation in order to arrive at the corresponding mass density, ρρ. In this paper we described our numerical methods to perform this reverse transform, and results are described for the specific qq from [3] including the computed ρρ, eigenvalues, and eigenfunctions for certain values of μ1μ1 and chosen initial conditions. Since we know what the eigenvalues are supposed to be for this qq, an analysis of the error in computing the eigenvalues is included. For comparison, similar calculations were done for other qq or ρρ in closed form.

## Acknowledgments

We would like to acknowledge the patience and guidance of Dr. Mark Embree and Dr. Steven Cox of Rice University. Many thanks are owed to Jeffrey Hokanson and our predecessors in the Physics of Strings PFUG, whose work made ours possible. This paper and the research herein conducted were made possible by grant funding from the National Science Foundation through Rice University's VIGRE program.

## Appendix

The codes used in this project are accessible through the featured link.

## References

1. Diacu, F. (2000). An Introduction to Differential Equations, Order and Chaos. W.H. Freeman and Company.
2. Hille, E. (1969). Lectures on Ordinary Differential Equations. Addison-Wesley Publishing Co.
3. Pöschel, J. and Trubowitz, E. (1986). Inverse Spectral Theory. Academic Press, Inc.
4. Trefethen, L. (2000). Spectral Methods in MATLAB. Society for Industrial and Applied Mathematics.
5. Truesdell, C. (1960). The Rational Mechanics of Flexible or Elastic Bodies. Introduction to Leonhardi Euleri Opera Omnia, Ser. II, Vol. X and XI. Leonhardi Euleri, Opera Omnia, Ser. II, Vol. XI (second part). Orell Füssli.

## Content actions

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

#### Collection to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

#### Module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks